import multi-element compatible pair style edip as edip/multi

This commit is contained in:
Axel Kohlmeyer 2017-05-16 17:40:04 -04:00
parent 06c151421c
commit 35e92733e9
18 changed files with 2146 additions and 8 deletions

View File

@ -1016,6 +1016,7 @@ package"_Section_start.html#start_3.
"dpd/fdt/energy"_pair_dpd_fdt.html,
"eam/cd (o)"_pair_eam.html,
"edip (o)"_pair_edip.html,
"edip/multi"_pair_edip.html,
"eff/cut"_pair_eff.html,
"exp6/rx"_pair_exp6_rx.html,
"gauss/cut"_pair_gauss.html,

View File

@ -12,6 +12,7 @@ pair_style edip command :h3
pair_style edip :pre
pair_style edip/omp :pre
pair_style edip/multi :pre
[Examples:]
@ -20,11 +21,14 @@ pair_coeff * * Si.edip Si
[Description:]
The {edip} style computes a 3-body "EDIP"_#EDIP potential which is
popular for modeling silicon materials where it can have advantages
over other models such as the "Stillinger-Weber"_pair_sw.html or
"Tersoff"_pair_tersoff.html potentials. In EDIP, the energy E of a
system of atoms is
The {edip} and {edip/multi} styles compute a 3-body "EDIP"_#EDIP
potential which is popular for modeling silicon materials where
it can have advantages over other models such as the
"Stillinger-Weber"_pair_sw.html or "Tersoff"_pair_tersoff.html
potentials. The {edip} style has been programmed for single element
potentials, while {edip/multi} supports multi-element EDIP runs.
In EDIP, the energy E of a system of atoms is
:c,image(Eqs/pair_edip.jpg)
@ -142,7 +146,7 @@ This pair style can only be used via the {pair} keyword of the
[Restrictions:]
This angle style can only be used if LAMMPS was built with the
This pair style can only be used if LAMMPS was built with the
USER-MISC package. See the "Making LAMMPS"_Section_start.html#start_3
section for more info on packages.
@ -151,7 +155,7 @@ for pair interactions.
The EDIP potential files provided with LAMMPS (see the potentials directory)
are parameterized for metal "units"_units.html.
You can use the SW potential with any LAMMPS units, but you would need
You can use the EDIP potential with any LAMMPS units, but you would need
to create your own EDIP potential file with coefficients listed in the
appropriate units if your simulation doesn't use "metal" units.

View File

@ -0,0 +1,26 @@
# DATE: 2011-09-15 CONTRIBUTOR: Unknown CITATION: Justo, Bazant, Kaxiras, Bulatov and Yip, Phys Rev B, 58, 2539 (1998)
# EDIP parameters for various elements and mixtures
# multiple entries can be added to this file, LAMMPS reads the ones it needs
# these entries are in LAMMPS "metal" units
# format of a single entry (one or more lines)
#
# element 1, element 2, element 3,
# A B cutoffA cutoffC alpha beta eta
# gamma lambda mu rho sigma Q0
# u1 u2 u3 u4
#
# units for each parameters:
# A , lambda are in eV
# B, cutoffA, cutoffC, gamma, sigma are in Angstrom
# alpha, beta, eta, mu, rho, Q0, u1-u4 are pure numbers
# Here are the original parameters in metal units, for Silicon from:
# J. F. Justo, M. Z. Bazant, E. Kaxiras, V. V. Bulatov, S. Yip
# Phys. Rev. B 58, 2539 (1998)
#
Si Si Si 7.9821730 1.5075463 3.1213820 2.5609104 3.1083847 0.0070975 0.2523244
1.1247945 1.4533108 0.6966326 1.2085196 0.5774108 312.1341346
-0.165799 32.557 0.286198 0.66

View File

@ -0,0 +1,38 @@
# DATE: 2017-05-16 CONTRIBUTOR: Laurent Pizzagalli CITATION: G. Lucas, M. Bertolus, and L. Pizzagalli, J. Phys. : Condens. Matter 22, 035802 (2010)
# element 1, element 2, element 3,
# A B cutoffA cutoffC alpha beta eta
# gamma lambda mu rho sigma Q0
# u1 u2 u3 u4
#
Si Si Si 5.488043 1.446435 2.941586 2.540193 3.066580 0.008593 0.589390
1.135256 2.417497 0.629131 1.343679 0.298443 208.924548
-0.165799 32.557 0.286198 0.66
C C C 10.222599 0.959814 2.212263 1.741598 1.962090 0.025661 0.275605
1.084183 3.633621 0.594236 2.827634 0.536561 289.305617
-0.165799 32.557 0.286198 0.66
C Si Si 7.535967 1.177019 2.534972 1.973974 2.507738 0.015347 0.432497
1.191567 3.025559 0.611684 2.061835 0.423863 249.115082
-0.165799 32.557000 0.286198 0.660000
Si C C 7.535967 1.177019 2.534972 1.973974 2.507738 0.015347 0.432497
1.191567 3.025559 0.611684 2.061835 0.423863 249.115082
-0.165799 32.557000 0.286198 0.660000
Si Si C 5.488043 1.446435 2.941586 2.540193 3.066580 0.008593 0.510944
1.135256 2.721528 0.620407 1.343679 0.298443 229.019815
-0.165799 32.557000 0.286198 0.660000
Si C Si 7.535967 1.177019 2.534972 1.973974 2.507738 0.015347 0.510944
1.191567 2.721528 0.620407 2.061835 0.423863 229.019815
-0.165799 32.557000 0.286198 0.660000
C C Si 10.222599 0.959814 2.212263 1.741598 1.962090 0.025661 0.354051
1.084183 3.329590 0.602960 2.827634 0.536561 269.210350
-0.165799 32.557000 0.286198 0.660000
C Si C 7.535967 1.177019 2.534972 1.973974 2.507738 0.015347 0.354051
1.191567 3.329590 0.602960 2.061835 0.423863 269.210350
-0.165799 32.557000 0.286198 0.660000

View File

@ -0,0 +1,138 @@
Position data for Silicon-Carbon system
128 atoms
2 atom types
-6.00 5.97232152 xlo xhi
-6.00 5.97232152 ylo yhi
-6.00 5.97232152 zlo zhi
Atoms
1 2 -2.9378454 -4.4592615 -4.8109196
2 2 5.6222143 -2.7335026 -1.7157569
3 2 -2.6614623 -5.5431059 1.6353686
4 2 -5.4326838 -4.6174577 5.9452279
5 2 5.8679239 -0.1120535 -3.5839373
6 2 -3.7174621 -0.6623311 -0.3714789
7 2 -5.0724728 -2.5671623 4.4103461
8 2 -3.3951436 0.9341126 4.9310702
9 2 -5.4347593 1.9523767 -5.6180938
10 2 -4.5884719 2.2904528 -1.0597739
11 2 -5.9058662 0.6212406 2.0127574
12 2 -4.7680660 0.1965740 4.3267764
13 2 -5.4228882 5.2569673 -4.5162920
14 2 -5.2683965 -5.9193658 -2.8648668
15 2 -2.8610884 1.0484664 2.0299077
16 2 -4.0711084 5.3133026 3.8009514
17 2 -0.1947147 -4.1677696 -5.6950931
18 2 -2.9892710 -3.1647368 -1.6173910
19 2 -0.9129311 -4.3819066 -0.1601859
20 2 -2.4513693 -5.2466501 4.8882912
21 2 -2.8879952 -0.1633446 -3.3401150
22 1 -4.6738762 -1.3807254 -2.2946777
23 2 -0.6973948 -1.4885343 0.6005156
24 1 -2.7392164 -2.4774843 0.2387186
25 2 -2.6551254 -2.7229952 2.6350264
26 1 -3.4644263 -4.6028144 3.3817786
27 2 0.7227614 -2.0709446 2.9214737
28 1 -2.1000577 -3.2131296 5.7273437
29 2 -3.1057649 2.3204819 -2.2725622
30 1 -2.2298751 0.7168389 -1.3107201
31 2 -1.8698261 1.4006751 0.7265108
32 1 -4.1103409 -0.7093340 1.9341753
33 2 -0.3505581 3.2707182 -0.2880656
34 1 -3.4045407 -1.4383961 4.3903527
35 2 -3.0940529 1.4132478 -5.3635505
36 1 -4.4560663 1.2072875 -3.7310176
37 2 -2.6061002 4.6373499 -4.6903941
38 1 -3.3477444 4.6768137 -2.6284678
39 2 0.8121697 4.8602418 -4.6710946
40 1 -2.5756922 3.3740738 -0.2136350
41 2 -0.3867976 5.8745611 -2.1119905
42 1 -1.6766249 1.3374292 3.8741477
43 2 -0.8770613 3.3735941 4.3846975
44 1 -1.8609254 3.3158245 -5.9786556
45 1 -5.2732321 -4.6073253 -0.9581754
46 1 -2.7888697 -5.6910152 -0.7922023
47 1 -2.4717165 4.5801880 2.5083210
48 1 -3.8819950 5.8456589 -5.7563384
49 2 2.2314782 -2.7729214 -5.2356862
50 2 0.2981976 -3.1385279 -3.1608167
51 2 2.8810785 -3.4658695 -0.5823196
52 2 0.2509625 -5.7595229 2.7389761
53 2 -0.2934120 -0.8029431 -3.3698507
54 1 -1.0075690 -2.0481922 -1.9419298
55 2 2.0729069 1.4922441 -2.3898096
56 1 1.1110944 -3.2004208 0.9491078
57 2 1.6774298 -0.7901860 2.5158773
58 1 -0.8342297 -4.3342518 2.0971458
59 2 3.2747406 -1.3107897 4.7884706
60 1 1.7126246 -3.3691471 4.5581012
61 2 0.4770605 1.7769008 -5.3339915
62 1 0.2944391 0.5892781 -2.2030106
63 2 2.2039275 3.1557557 -2.0276796
64 1 -0.0404494 0.4767818 1.0396418
65 2 1.1395867 2.3763443 2.3481007
66 1 -0.9738374 -1.6325161 3.7538567
67 2 -0.3291998 0.2996990 5.2770809
68 1 -1.6185604 -0.3964274 -5.1771220
69 2 2.5999949 -5.1977715 5.8230717
70 1 -1.6270675 2.3210900 -3.6299941
71 2 3.6532700 4.9282597 -5.4319276
72 1 0.0788934 4.0241037 -2.5011530
73 2 2.8556507 2.6168653 2.1125546
74 1 0.9738989 2.6255364 4.3412121
75 2 3.7452938 3.4521356 4.5946426
76 1 2.0805182 4.7039015 5.3280260
77 1 -1.0324174 -5.8155041 -4.3265820
78 1 0.7622442 -4.3631629 -1.3156572
79 1 0.3263684 3.9937357 1.6172321
80 1 -0.4350105 -5.7997058 4.5959134
81 2 3.9161132 -4.6052788 -3.3191717
82 2 1.9240657 5.7345079 -1.9754251
83 2 -5.9794488 -4.2369359 1.8646522
84 2 4.3339975 -4.4845227 5.3737440
85 2 2.2755456 -0.6327737 -5.7931837
86 1 1.8728190 -1.5504906 -3.4560010
87 2 3.4558100 -1.1054068 -1.8333071
88 1 4.3788172 -1.9466494 -0.3284637
89 2 2.5999235 -3.7548996 2.5740569
90 1 3.9983910 -4.4856603 1.1968663
91 2 -5.7295580 -2.1475672 -5.9963645
92 1 4.2664051 -2.6988975 -5.8005478
93 2 4.5254685 2.2906832 -3.4765798
94 1 2.3603088 1.3416442 -4.4173836
95 2 4.7767057 1.4061217 -0.7524620
96 1 1.8072666 -0.7835973 -0.4581995
97 2 4.4745018 0.3736224 2.1068274
98 1 3.6081170 -1.7315713 2.4019053
99 2 4.6281423 -0.2865409 4.4756524
100 1 1.7975239 0.2893530 4.2330830
101 2 5.8341452 4.4986472 -5.9664541
102 1 3.2401308 4.1655227 -3.5070029
103 2 4.8720339 4.8709982 -2.3364366
104 1 3.5526476 1.2262752 0.6926826
105 2 -5.8173342 4.5420479 1.5578881
106 1 3.9683224 1.5441137 3.8284375
107 2 -5.5349308 1.9067049 3.7504113
108 1 4.4728615 2.6415574 -5.5952809
109 1 1.7000950 -4.8115440 -4.1953920
110 1 1.7221527 4.1878404 -0.3712681
111 1 3.9218156 4.5935583 1.3263407
112 1 3.1310195 -5.8922481 3.6001155
113 1 4.7558719 -2.2877771 -3.4742052
114 1 -5.5050300 -2.7027381 0.8748867
115 1 5.8418594 -4.6064370 3.8714113
116 1 -4.7516868 -3.1691984 -4.4099768
117 1 3.9404971 0.7188702 -2.2898786
118 1 -5.6869740 0.2042380 -0.1916738
119 1 5.8949589 -1.2422560 3.1201292
120 1 5.9675804 -0.0712572 5.8964022
121 1 -5.6208517 3.3600036 -2.9493510
122 1 5.2065263 3.4517912 -0.3800894
123 1 -4.6994522 2.5489583 1.8297431
124 1 -4.0758407 3.0726196 5.0647973
125 1 4.1587591 -5.0896820 -1.1443498
126 1 -4.6963753 -5.7429833 1.1357818
127 1 5.5994192 4.6887008 3.5948264
128 1 5.0988369 -5.3774409 -4.9051267

View File

@ -0,0 +1,72 @@
units metal
atom_style atomic
atom_modify map array
boundary p p p
atom_modify sort 0 0.0
# temperature
variable t equal 1800.0
# coordination number cutoff
variable r equal 2.835
# minimization parameters
variable etol equal 1.0e-5
variable ftol equal 1.0e-5
variable maxiter equal 100
variable maxeval equal 100
variable dmax equal 1.0e-1
# diamond unit cell
variable a equal 5.431
lattice custom $a &
a1 1.0 0.0 0.0 &
a2 0.0 1.0 0.0 &
a3 0.0 0.0 1.0 &
basis 0.0 0.0 0.0 &
basis 0.0 0.5 0.5 &
basis 0.5 0.0 0.5 &
basis 0.5 0.5 0.0 &
basis 0.25 0.25 0.25 &
basis 0.25 0.75 0.75 &
basis 0.75 0.25 0.75 &
basis 0.75 0.75 0.25
region myreg block 0 4 &
0 4 &
0 4
create_box 1 myreg
create_atoms 1 region myreg
mass 1 28.06
group Si type 1
velocity all create $t 5287287 mom yes rot yes dist gaussian
# make a vacancy
group del id 300
delete_atoms group del
pair_style edip
pair_coeff * * Si.edip Si
thermo 10
fix 1 all nvt temp $t $t 0.1
timestep 1.0e-3
neighbor 1.0 bin
neigh_modify every 1 delay 10 check yes
# equilibrate
run 500

View File

@ -0,0 +1,72 @@
units metal
atom_style atomic
atom_modify map array
boundary p p p
atom_modify sort 0 0.0
# temperature
variable t equal 1800.0
# coordination number cutoff
variable r equal 2.835
# minimization parameters
variable etol equal 1.0e-5
variable ftol equal 1.0e-5
variable maxiter equal 100
variable maxeval equal 100
variable dmax equal 1.0e-1
# diamond unit cell
variable a equal 5.431
lattice custom $a &
a1 1.0 0.0 0.0 &
a2 0.0 1.0 0.0 &
a3 0.0 0.0 1.0 &
basis 0.0 0.0 0.0 &
basis 0.0 0.5 0.5 &
basis 0.5 0.0 0.5 &
basis 0.5 0.5 0.0 &
basis 0.25 0.25 0.25 &
basis 0.25 0.75 0.75 &
basis 0.75 0.25 0.75 &
basis 0.75 0.75 0.25
region myreg block 0 4 &
0 4 &
0 4
create_box 1 myreg
create_atoms 1 region myreg
mass 1 28.06
group Si type 1
velocity all create $t 5287287 mom yes rot yes dist gaussian
# make a vacancy
group del id 300
delete_atoms group del
pair_style edip/multi
pair_coeff * * Si.edip Si
thermo 10
fix 1 all nvt temp $t $t 0.1
timestep 1.0e-3
neighbor 1.0 bin
neigh_modify every 1 delay 10 check yes
# equilibrate
run 500

View File

@ -0,0 +1,33 @@
# Test of MEAM potential for SiC system
units metal
boundary p p p
atom_style atomic
read_data data.SiC
pair_style edip/multi
pair_coeff * * SiC.edip Si C
mass 1 28.085
mass 2 12.001
neighbor 1.0 bin
neigh_modify delay 1
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

View File

@ -0,0 +1,167 @@
LAMMPS (4 May 2017)
using 1 OpenMP thread(s) per MPI task
units metal
atom_style atomic
atom_modify map array
boundary p p p
atom_modify sort 0 0.0
# temperature
variable t equal 1800.0
# coordination number cutoff
variable r equal 2.835
# minimization parameters
variable etol equal 1.0e-5
variable ftol equal 1.0e-5
variable maxiter equal 100
variable maxeval equal 100
variable dmax equal 1.0e-1
# diamond unit cell
variable a equal 5.431
lattice custom $a a1 1.0 0.0 0.0 a2 0.0 1.0 0.0 a3 0.0 0.0 1.0 basis 0.0 0.0 0.0 basis 0.0 0.5 0.5 basis 0.5 0.0 0.5 basis 0.5 0.5 0.0 basis 0.25 0.25 0.25 basis 0.25 0.75 0.75 basis 0.75 0.25 0.75 basis 0.75 0.75 0.25
lattice custom 5.431 a1 1.0 0.0 0.0 a2 0.0 1.0 0.0 a3 0.0 0.0 1.0 basis 0.0 0.0 0.0 basis 0.0 0.5 0.5 basis 0.5 0.0 0.5 basis 0.5 0.5 0.0 basis 0.25 0.25 0.25 basis 0.25 0.75 0.75 basis 0.75 0.25 0.75 basis 0.75 0.75 0.25
Lattice spacing in x,y,z = 5.431 5.431 5.431
region myreg block 0 4 0 4 0 4
create_box 1 myreg
Created orthogonal box = (0 0 0) to (21.724 21.724 21.724)
1 by 1 by 1 MPI processor grid
create_atoms 1 region myreg
Created 512 atoms
mass 1 28.06
group Si type 1
512 atoms in group Si
velocity all create $t 5287287 mom yes rot yes dist gaussian
velocity all create 1800 5287287 mom yes rot yes dist gaussian
# make a vacancy
group del id 300
1 atoms in group del
delete_atoms group del
Deleted 1 atoms, new total = 511
pair_style edip/multi
pair_coeff * * Si.edip Si
Reading potential file Si.edip with DATE: 2011-09-15
thermo 10
fix 1 all nvt temp $t $t 0.1
fix 1 all nvt temp 1800 $t 0.1
fix 1 all nvt temp 1800 1800 0.1
timestep 1.0e-3
neighbor 1.0 bin
neigh_modify every 1 delay 10 check yes
# equilibrate
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 = 4.12138
ghost atom cutoff = 4.12138
binsize = 2.06069, bins = 11 11 11
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair edip/multi, perpetual
attributes: full, newton on
pair build: full/bin/atomonly
stencil: full/bin/3d
bin: standard
Per MPI rank memory allocation (min/avg/max) = 2.979 | 2.979 | 2.979 Mbytes
Step Temp E_pair E_mol TotEng Press
0 1802.5039 -2372.6618 0 -2253.8359 12261.807
10 952.62744 -2316.428 0 -2253.6283 723.08194
20 549.13801 -2289.442 0 -2253.2413 -2444.5204
30 1047.0106 -2321.1523 0 -2252.1305 9013.201
40 663.46141 -2294.2083 0 -2250.4711 2942.5348
50 504.74535 -2282.849 0 -2249.5748 -461.44909
60 1019.2173 -2315.5639 0 -2248.3744 7706.4286
70 844.51195 -2302.5251 0 -2246.8526 3116.8302
80 814.90407 -2299.3372 0 -2245.6166 794.77455
90 1269.5636 -2327.4775 0 -2243.7845 7729.3968
100 977.61563 -2306.1118 0 -2241.6647 2969.9939
110 843.08539 -2295.6547 0 -2240.0763 1393.4039
120 1161.6968 -2314.6587 0 -2238.0766 7398.3492
130 918.19451 -2296.4321 0 -2235.9022 2537.3997
140 881.42548 -2292.2768 0 -2234.1709 1550.3339
150 1231.1005 -2313.1054 0 -2231.9479 8112.7566
160 967.01862 -2293.332 0 -2229.5836 3422.9627
170 833.51248 -2282.7489 0 -2227.8015 43.991459
180 1240.8488 -2307.3633 0 -2225.5632 6557.8651
190 1126.4621 -2297.1922 0 -2222.9328 4289.0067
200 947.59571 -2283.29 0 -2220.822 586.2811
210 1228.153 -2299.4702 0 -2218.5071 5315.0425
220 1215.4104 -2295.9408 0 -2215.8176 4870.3417
230 1112.436 -2286.7552 0 -2213.4204 2527.1879
240 1300.081 -2296.6013 0 -2210.8965 5738.3708
250 1192.5738 -2286.8463 0 -2208.2286 4076.49
260 1004.7055 -2272.1753 0 -2205.9424 359.37589
270 1241.2018 -2285.3632 0 -2203.5399 4160.5763
280 1360.1974 -2290.325 0 -2200.6572 5802.3902
290 1151.9365 -2273.9467 0 -2198.008 1418.8887
300 1174.3518 -2273.0089 0 -2195.5925 1998.229
310 1329.2727 -2280.5049 0 -2192.8757 4721.7297
320 1284.4414 -2274.7519 0 -2190.0781 2985.4674
330 1328.3761 -2274.9545 0 -2187.3844 4543.2109
340 1446.3847 -2279.8693 0 -2184.5198 6254.4059
350 1366.2165 -2271.7475 0 -2181.6828 3637.8335
360 1358.9609 -2268.5982 0 -2179.0118 3049.5798
370 1552.208 -2278.4802 0 -2176.1545 6334.0058
380 1562.5295 -2276.1793 0 -2173.1732 5787.5547
390 1415.5498 -2263.7824 0 -2170.4655 3438.5766
400 1323.1568 -2255.1641 0 -2167.938 2427.2294
410 1260.7186 -2248.5373 0 -2165.4273 1208.6299
420 1282.1118 -2247.3718 0 -2162.8516 462.65374
430 1451.944 -2255.7551 0 -2160.0391 2037.8025
440 1568.9415 -2260.417 0 -2156.9882 3531.1602
450 1565.8262 -2257.2396 0 -2154.0162 2586.7886
460 1677.7143 -2261.7214 0 -2151.122 4112.9756
470 1762.9071 -2264.4244 0 -2148.2089 5053.2139
480 1704.5898 -2257.8678 0 -2145.4967 4077.4626
490 1731.2619 -2257.1048 0 -2142.9753 4710.5263
500 1723.9777 -2254.161 0 -2140.5118 4760.7295
Loop time of 0.679564 on 1 procs for 500 steps with 511 atoms
Performance: 63.570 ns/day, 0.378 hours/ns, 735.765 timesteps/s
99.7% CPU use with 1 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 0.65181 | 0.65181 | 0.65181 | 0.0 | 95.92
Neigh | 0.013857 | 0.013857 | 0.013857 | 0.0 | 2.04
Comm | 0.0033884 | 0.0033884 | 0.0033884 | 0.0 | 0.50
Output | 0.00070739 | 0.00070739 | 0.00070739 | 0.0 | 0.10
Modify | 0.0083694 | 0.0083694 | 0.0083694 | 0.0 | 1.23
Other | | 0.001432 | | | 0.21
Nlocal: 511 ave 511 max 511 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 845 ave 845 max 845 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 0 ave 0 max 0 min
Histogram: 1 0 0 0 0 0 0 0 0 0
FullNghs: 7902 ave 7902 max 7902 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 7902
Ave neighs/atom = 15.4638
Neighbor list builds = 19
Dangerous builds = 0
Total wall time: 0:00:00

View File

@ -0,0 +1,167 @@
LAMMPS (4 May 2017)
using 1 OpenMP thread(s) per MPI task
units metal
atom_style atomic
atom_modify map array
boundary p p p
atom_modify sort 0 0.0
# temperature
variable t equal 1800.0
# coordination number cutoff
variable r equal 2.835
# minimization parameters
variable etol equal 1.0e-5
variable ftol equal 1.0e-5
variable maxiter equal 100
variable maxeval equal 100
variable dmax equal 1.0e-1
# diamond unit cell
variable a equal 5.431
lattice custom $a a1 1.0 0.0 0.0 a2 0.0 1.0 0.0 a3 0.0 0.0 1.0 basis 0.0 0.0 0.0 basis 0.0 0.5 0.5 basis 0.5 0.0 0.5 basis 0.5 0.5 0.0 basis 0.25 0.25 0.25 basis 0.25 0.75 0.75 basis 0.75 0.25 0.75 basis 0.75 0.75 0.25
lattice custom 5.431 a1 1.0 0.0 0.0 a2 0.0 1.0 0.0 a3 0.0 0.0 1.0 basis 0.0 0.0 0.0 basis 0.0 0.5 0.5 basis 0.5 0.0 0.5 basis 0.5 0.5 0.0 basis 0.25 0.25 0.25 basis 0.25 0.75 0.75 basis 0.75 0.25 0.75 basis 0.75 0.75 0.25
Lattice spacing in x,y,z = 5.431 5.431 5.431
region myreg block 0 4 0 4 0 4
create_box 1 myreg
Created orthogonal box = (0 0 0) to (21.724 21.724 21.724)
1 by 2 by 2 MPI processor grid
create_atoms 1 region myreg
Created 512 atoms
mass 1 28.06
group Si type 1
512 atoms in group Si
velocity all create $t 5287287 mom yes rot yes dist gaussian
velocity all create 1800 5287287 mom yes rot yes dist gaussian
# make a vacancy
group del id 300
1 atoms in group del
delete_atoms group del
Deleted 1 atoms, new total = 511
pair_style edip/multi
pair_coeff * * Si.edip Si
Reading potential file Si.edip with DATE: 2011-09-15
thermo 10
fix 1 all nvt temp $t $t 0.1
fix 1 all nvt temp 1800 $t 0.1
fix 1 all nvt temp 1800 1800 0.1
timestep 1.0e-3
neighbor 1.0 bin
neigh_modify every 1 delay 10 check yes
# equilibrate
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 = 4.12138
ghost atom cutoff = 4.12138
binsize = 2.06069, bins = 11 11 11
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair edip/multi, perpetual
attributes: full, newton on
pair build: full/bin/atomonly
stencil: full/bin/3d
bin: standard
Per MPI rank memory allocation (min/avg/max) = 2.955 | 2.955 | 2.955 Mbytes
Step Temp E_pair E_mol TotEng Press
0 1802.3816 -2372.6618 0 -2253.844 12260.967
10 938.75954 -2315.5185 0 -2253.6329 558.21646
20 534.27233 -2288.4721 0 -2253.2514 -2710.768
30 1043.7796 -2320.9485 0 -2252.1398 8679.4381
40 658.0916 -2293.8597 0 -2250.4765 2165.3742
50 517.93009 -2283.7238 0 -2249.5805 -1124.9373
60 1063.3594 -2318.4409 0 -2248.3414 7277.8526
70 868.14006 -2304.0134 0 -2246.7832 2050.2848
80 826.37805 -2300.0187 0 -2245.5416 91.099408
90 1289.6772 -2328.7151 0 -2243.6961 8180.7423
100 976.36208 -2305.9371 0 -2241.5727 3614.0499
110 810.81713 -2293.4705 0 -2240.0193 1359.368
120 1165.707 -2314.9026 0 -2238.056 7336.45
130 929.81245 -2297.139 0 -2235.8432 2793.8451
140 804.47874 -2287.2074 0 -2234.174 704.92455
150 1182.4141 -2310.0266 0 -2232.0787 7822.2337
160 979.92391 -2294.2969 0 -2229.6977 3206.7458
170 830.14748 -2282.6079 0 -2227.8824 -296.87377
180 1271.1133 -2309.4274 0 -2225.6322 7199.614
190 1209.6006 -2302.6407 0 -2222.9006 5528.3784
200 954.67693 -2283.6621 0 -2220.7273 47.02795
210 1260.814 -2301.5582 0 -2218.442 4829.788
220 1274.9954 -2299.7285 0 -2215.6774 5518.0597
230 1048.0074 -2282.398 0 -2213.3106 1754.4144
240 1261.7072 -2294.1108 0 -2210.9356 5233.2712
250 1272.6178 -2292.0793 0 -2208.1849 4795.9325
260 989.14205 -2271.0278 0 -2205.8209 -820.1828
270 1212.0445 -2283.4212 0 -2203.52 3395.8634
280 1391.9572 -2292.3809 0 -2200.6194 6666.2451
290 1093.1204 -2270.0421 0 -2197.9807 206.94523
300 1159.4831 -2272.102 0 -2195.6657 778.53806
310 1407.3528 -2285.6228 0 -2192.8463 5223.048
320 1236.7163 -2271.5389 0 -2190.0113 1865.3943
330 1258.8275 -2270.4611 0 -2187.4758 2333.3209
340 1507.9519 -2283.9906 0 -2184.5824 6775.5456
350 1366.5116 -2271.7287 0 -2181.6446 3432.115
360 1305.2829 -2265.1092 0 -2179.0614 1498.4073
370 1581.4335 -2280.4645 0 -2176.2122 6518.5597
380 1589.5319 -2277.9428 0 -2173.1567 6334.6506
390 1402.6781 -2262.9323 0 -2170.464 3278.3038
400 1374.9587 -2258.5717 0 -2167.9307 3608.7284
410 1295.7416 -2250.7752 0 -2165.3565 1877.5222
420 1278.6727 -2247.1099 0 -2162.8164 1599.4181
430 1508.1328 -2259.4245 0 -2160.0044 4300.2224
440 1624.2957 -2263.9806 0 -2156.9026 4432.625
450 1597.3356 -2259.263 0 -2153.9624 3370.3816
460 1772.0922 -2267.9106 0 -2151.0895 5788.3214
470 1806.4047 -2267.304 0 -2148.221 5950.1166
480 1593.0406 -2250.7469 0 -2145.7294 2518.0576
490 1660.9767 -2252.894 0 -2143.398 4282.1643
500 1714.283 -2253.9295 0 -2140.9194 5740.0247
Loop time of 0.205398 on 4 procs for 500 steps with 511 atoms
Performance: 210.324 ns/day, 0.114 hours/ns, 2434.304 timesteps/s
99.0% 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.16285 | 0.1688 | 0.17446 | 1.1 | 82.18
Neigh | 0.0035172 | 0.0036234 | 0.0038214 | 0.2 | 1.76
Comm | 0.018727 | 0.024851 | 0.030996 | 2.9 | 12.10
Output | 0.0013061 | 0.0014012 | 0.0015635 | 0.3 | 0.68
Modify | 0.0046582 | 0.0048603 | 0.0050988 | 0.2 | 2.37
Other | | 0.001861 | | | 0.91
Nlocal: 127.75 ave 131 max 124 min
Histogram: 1 0 1 0 0 0 0 0 1 1
Nghost: 433.75 ave 441 max 426 min
Histogram: 1 0 1 0 0 0 0 0 1 1
Neighs: 0 ave 0 max 0 min
Histogram: 4 0 0 0 0 0 0 0 0 0
FullNghs: 1979.5 ave 2040 max 1895 min
Histogram: 1 0 0 0 1 0 0 0 0 2
Total # of neighbors = 7918
Ave neighs/atom = 15.4951
Neighbor list builds = 19
Dangerous builds = 0
Total wall time: 0:00:00

View File

@ -0,0 +1,167 @@
LAMMPS (4 May 2017)
using 1 OpenMP thread(s) per MPI task
units metal
atom_style atomic
atom_modify map array
boundary p p p
atom_modify sort 0 0.0
# temperature
variable t equal 1800.0
# coordination number cutoff
variable r equal 2.835
# minimization parameters
variable etol equal 1.0e-5
variable ftol equal 1.0e-5
variable maxiter equal 100
variable maxeval equal 100
variable dmax equal 1.0e-1
# diamond unit cell
variable a equal 5.431
lattice custom $a a1 1.0 0.0 0.0 a2 0.0 1.0 0.0 a3 0.0 0.0 1.0 basis 0.0 0.0 0.0 basis 0.0 0.5 0.5 basis 0.5 0.0 0.5 basis 0.5 0.5 0.0 basis 0.25 0.25 0.25 basis 0.25 0.75 0.75 basis 0.75 0.25 0.75 basis 0.75 0.75 0.25
lattice custom 5.431 a1 1.0 0.0 0.0 a2 0.0 1.0 0.0 a3 0.0 0.0 1.0 basis 0.0 0.0 0.0 basis 0.0 0.5 0.5 basis 0.5 0.0 0.5 basis 0.5 0.5 0.0 basis 0.25 0.25 0.25 basis 0.25 0.75 0.75 basis 0.75 0.25 0.75 basis 0.75 0.75 0.25
Lattice spacing in x,y,z = 5.431 5.431 5.431
region myreg block 0 4 0 4 0 4
create_box 1 myreg
Created orthogonal box = (0 0 0) to (21.724 21.724 21.724)
1 by 1 by 1 MPI processor grid
create_atoms 1 region myreg
Created 512 atoms
mass 1 28.06
group Si type 1
512 atoms in group Si
velocity all create $t 5287287 mom yes rot yes dist gaussian
velocity all create 1800 5287287 mom yes rot yes dist gaussian
# make a vacancy
group del id 300
1 atoms in group del
delete_atoms group del
Deleted 1 atoms, new total = 511
pair_style edip
pair_coeff * * Si.edip Si
Reading potential file Si.edip with DATE: 2011-09-15
thermo 10
fix 1 all nvt temp $t $t 0.1
fix 1 all nvt temp 1800 $t 0.1
fix 1 all nvt temp 1800 1800 0.1
timestep 1.0e-3
neighbor 1.0 bin
neigh_modify every 1 delay 10 check yes
# equilibrate
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 = 4.12138
ghost atom cutoff = 4.12138
binsize = 2.06069, bins = 11 11 11
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair edip, perpetual
attributes: full, newton on
pair build: full/bin/atomonly
stencil: full/bin/3d
bin: standard
Per MPI rank memory allocation (min/avg/max) = 2.979 | 2.979 | 2.979 Mbytes
Step Temp E_pair E_mol TotEng Press
0 1802.5039 -2372.6618 0 -2253.8359 12261.807
10 952.62744 -2316.428 0 -2253.6283 723.08283
20 549.138 -2289.442 0 -2253.2413 -2444.5194
30 1047.0106 -2321.1522 0 -2252.1305 9013.2015
40 663.46143 -2294.2083 0 -2250.4711 2942.5358
50 504.74533 -2282.849 0 -2249.5748 -461.44817
60 1019.2173 -2315.5639 0 -2248.3744 7706.429
70 844.51197 -2302.5251 0 -2246.8526 3116.8313
80 814.90406 -2299.3372 0 -2245.6165 794.77536
90 1269.5635 -2327.4775 0 -2243.7845 7729.3971
100 977.61566 -2306.1118 0 -2241.6647 2969.9952
110 843.08538 -2295.6547 0 -2240.0763 1393.4046
120 1161.6968 -2314.6587 0 -2238.0766 7398.3495
130 918.19453 -2296.4321 0 -2235.9022 2537.4011
140 881.42546 -2292.2768 0 -2234.1709 1550.3345
150 1231.1005 -2313.1054 0 -2231.9479 8112.7568
160 967.01865 -2293.332 0 -2229.5836 3422.964
170 833.51246 -2282.7489 0 -2227.8015 43.99251
180 1240.8487 -2307.3633 0 -2225.5632 6557.8652
190 1126.4621 -2297.1922 0 -2222.9328 4289.0083
200 947.5957 -2283.29 0 -2220.8219 586.28203
210 1228.153 -2299.4702 0 -2218.5071 5315.0427
220 1215.4104 -2295.9407 0 -2215.8176 4870.343
230 1112.436 -2286.7552 0 -2213.4204 2527.1887
240 1300.081 -2296.6013 0 -2210.8965 5738.3711
250 1192.5739 -2286.8463 0 -2208.2286 4076.4913
260 1004.7055 -2272.1753 0 -2205.9424 359.3769
270 1241.2018 -2285.3632 0 -2203.5399 4160.5764
280 1360.1974 -2290.325 0 -2200.6572 5802.3912
290 1151.9366 -2273.9467 0 -2198.008 1418.8905
300 1174.3518 -2273.0089 0 -2195.5925 1998.2297
310 1329.2726 -2280.5049 0 -2192.8757 4721.7304
320 1284.4414 -2274.7519 0 -2190.0781 2985.4687
330 1328.3761 -2274.9545 0 -2187.3844 4543.2115
340 1446.3847 -2279.8693 0 -2184.5198 6254.4071
350 1366.2165 -2271.7475 0 -2181.6828 3637.8351
360 1358.9609 -2268.5982 0 -2179.0118 3049.5811
370 1552.2079 -2278.4802 0 -2176.1545 6334.0061
380 1562.5295 -2276.1793 0 -2173.1731 5787.5565
390 1415.5498 -2263.7823 0 -2170.4655 3438.5782
400 1323.1568 -2255.1641 0 -2167.938 2427.2311
410 1260.7186 -2248.5373 0 -2165.4273 1208.6316
420 1282.1118 -2247.3718 0 -2162.8516 462.65508
430 1451.9439 -2255.7551 0 -2160.0391 2037.8027
440 1568.9415 -2260.417 0 -2156.9882 3531.1613
450 1565.8261 -2257.2396 0 -2154.0161 2586.7896
460 1677.7143 -2261.7214 0 -2151.122 4112.976
470 1762.9071 -2264.4244 0 -2148.2089 5053.2148
480 1704.5898 -2257.8678 0 -2145.4966 4077.4649
490 1731.2619 -2257.1048 0 -2142.9753 4710.5276
500 1723.9777 -2254.161 0 -2140.5118 4760.7316
Loop time of 0.312472 on 1 procs for 500 steps with 511 atoms
Performance: 138.252 ns/day, 0.174 hours/ns, 1600.143 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.28525 | 0.28525 | 0.28525 | 0.0 | 91.29
Neigh | 0.013753 | 0.013753 | 0.013753 | 0.0 | 4.40
Comm | 0.0033333 | 0.0033333 | 0.0033333 | 0.0 | 1.07
Output | 0.00071096 | 0.00071096 | 0.00071096 | 0.0 | 0.23
Modify | 0.008044 | 0.008044 | 0.008044 | 0.0 | 2.57
Other | | 0.001385 | | | 0.44
Nlocal: 511 ave 511 max 511 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 845 ave 845 max 845 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 0 ave 0 max 0 min
Histogram: 1 0 0 0 0 0 0 0 0 0
FullNghs: 7902 ave 7902 max 7902 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 7902
Ave neighs/atom = 15.4638
Neighbor list builds = 19
Dangerous builds = 0
Total wall time: 0:00:00

View File

@ -0,0 +1,167 @@
LAMMPS (4 May 2017)
using 1 OpenMP thread(s) per MPI task
units metal
atom_style atomic
atom_modify map array
boundary p p p
atom_modify sort 0 0.0
# temperature
variable t equal 1800.0
# coordination number cutoff
variable r equal 2.835
# minimization parameters
variable etol equal 1.0e-5
variable ftol equal 1.0e-5
variable maxiter equal 100
variable maxeval equal 100
variable dmax equal 1.0e-1
# diamond unit cell
variable a equal 5.431
lattice custom $a a1 1.0 0.0 0.0 a2 0.0 1.0 0.0 a3 0.0 0.0 1.0 basis 0.0 0.0 0.0 basis 0.0 0.5 0.5 basis 0.5 0.0 0.5 basis 0.5 0.5 0.0 basis 0.25 0.25 0.25 basis 0.25 0.75 0.75 basis 0.75 0.25 0.75 basis 0.75 0.75 0.25
lattice custom 5.431 a1 1.0 0.0 0.0 a2 0.0 1.0 0.0 a3 0.0 0.0 1.0 basis 0.0 0.0 0.0 basis 0.0 0.5 0.5 basis 0.5 0.0 0.5 basis 0.5 0.5 0.0 basis 0.25 0.25 0.25 basis 0.25 0.75 0.75 basis 0.75 0.25 0.75 basis 0.75 0.75 0.25
Lattice spacing in x,y,z = 5.431 5.431 5.431
region myreg block 0 4 0 4 0 4
create_box 1 myreg
Created orthogonal box = (0 0 0) to (21.724 21.724 21.724)
1 by 2 by 2 MPI processor grid
create_atoms 1 region myreg
Created 512 atoms
mass 1 28.06
group Si type 1
512 atoms in group Si
velocity all create $t 5287287 mom yes rot yes dist gaussian
velocity all create 1800 5287287 mom yes rot yes dist gaussian
# make a vacancy
group del id 300
1 atoms in group del
delete_atoms group del
Deleted 1 atoms, new total = 511
pair_style edip
pair_coeff * * Si.edip Si
Reading potential file Si.edip with DATE: 2011-09-15
thermo 10
fix 1 all nvt temp $t $t 0.1
fix 1 all nvt temp 1800 $t 0.1
fix 1 all nvt temp 1800 1800 0.1
timestep 1.0e-3
neighbor 1.0 bin
neigh_modify every 1 delay 10 check yes
# equilibrate
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 = 4.12138
ghost atom cutoff = 4.12138
binsize = 2.06069, bins = 11 11 11
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair edip, perpetual
attributes: full, newton on
pair build: full/bin/atomonly
stencil: full/bin/3d
bin: standard
Per MPI rank memory allocation (min/avg/max) = 2.955 | 2.955 | 2.955 Mbytes
Step Temp E_pair E_mol TotEng Press
0 1802.3816 -2372.6618 0 -2253.8439 12260.967
10 938.75954 -2315.5185 0 -2253.6329 558.21736
20 534.27232 -2288.4721 0 -2253.2514 -2710.767
30 1043.7796 -2320.9485 0 -2252.1398 8679.4385
40 658.09162 -2293.8597 0 -2250.4765 2165.3752
50 517.93008 -2283.7238 0 -2249.5805 -1124.9362
60 1063.3594 -2318.4409 0 -2248.3414 7277.853
70 868.14007 -2304.0133 0 -2246.7832 2050.2859
80 826.37803 -2300.0187 0 -2245.5416 91.100098
90 1289.6772 -2328.7151 0 -2243.6961 8180.7427
100 976.36211 -2305.9371 0 -2241.5727 3614.0511
110 810.81711 -2293.4705 0 -2240.0193 1359.3687
120 1165.707 -2314.9026 0 -2238.056 7336.4505
130 929.81248 -2297.139 0 -2235.8432 2793.8463
140 804.47872 -2287.2074 0 -2234.174 704.92524
150 1182.414 -2310.0266 0 -2232.0787 7822.2339
160 979.92395 -2294.2969 0 -2229.6977 3206.7474
170 830.14746 -2282.6079 0 -2227.8824 -296.87288
180 1271.1133 -2309.4274 0 -2225.6322 7199.614
190 1209.6006 -2302.6407 0 -2222.9006 5528.3799
200 954.67692 -2283.6621 0 -2220.7272 47.02925
210 1260.814 -2301.5582 0 -2218.442 4829.7879
220 1274.9954 -2299.7285 0 -2215.6774 5518.0611
230 1048.0074 -2282.398 0 -2213.3106 1754.4157
240 1261.7071 -2294.1107 0 -2210.9356 5233.2714
250 1272.6179 -2292.0793 0 -2208.1849 4795.934
260 989.14207 -2271.0278 0 -2205.8209 -820.18098
270 1212.0444 -2283.4212 0 -2203.52 3395.8631
280 1391.9572 -2292.3809 0 -2200.6194 6666.2464
290 1093.1205 -2270.0421 0 -2197.9807 206.94752
300 1159.483 -2272.102 0 -2195.6657 778.53823
310 1407.3528 -2285.6227 0 -2192.8463 5223.0487
320 1236.7164 -2271.5389 0 -2190.0112 1865.3963
330 1258.8275 -2270.4611 0 -2187.4758 2333.321
340 1507.9519 -2283.9906 0 -2184.5824 6775.546
350 1366.5116 -2271.7287 0 -2181.6446 3432.1175
360 1305.2828 -2265.1091 0 -2179.0614 1498.4079
370 1581.4334 -2280.4645 0 -2176.2122 6518.5598
380 1589.5319 -2277.9428 0 -2173.1566 6334.6527
390 1402.6782 -2262.9323 0 -2170.464 3278.3048
400 1374.9587 -2258.5717 0 -2167.9307 3608.7293
410 1295.7416 -2250.7752 0 -2165.3565 1877.5245
420 1278.6727 -2247.1099 0 -2162.8164 1599.4189
430 1508.1328 -2259.4245 0 -2160.0044 4300.2235
440 1624.2957 -2263.9806 0 -2156.9026 4432.6267
450 1597.3356 -2259.263 0 -2153.9623 3370.3829
460 1772.0921 -2267.9105 0 -2151.0895 5788.3219
470 1806.4047 -2267.304 0 -2148.221 5950.1188
480 1593.0406 -2250.7469 0 -2145.7294 2518.0601
490 1660.9766 -2252.894 0 -2143.398 4282.1654
500 1714.2831 -2253.9295 0 -2140.9194 5740.0268
Loop time of 0.109584 on 4 procs for 500 steps with 511 atoms
Performance: 394.220 ns/day, 0.061 hours/ns, 4562.726 timesteps/s
99.0% 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.074678 | 0.077817 | 0.084705 | 1.4 | 71.01
Neigh | 0.0036662 | 0.0037943 | 0.0039661 | 0.2 | 3.46
Comm | 0.013665 | 0.020312 | 0.023178 | 2.7 | 18.54
Output | 0.0010247 | 0.0010931 | 0.0012922 | 0.3 | 1.00
Modify | 0.0043213 | 0.0047521 | 0.0051889 | 0.6 | 4.34
Other | | 0.001814 | | | 1.66
Nlocal: 127.75 ave 131 max 124 min
Histogram: 1 0 1 0 0 0 0 0 1 1
Nghost: 433.75 ave 441 max 426 min
Histogram: 1 0 1 0 0 0 0 0 1 1
Neighs: 0 ave 0 max 0 min
Histogram: 4 0 0 0 0 0 0 0 0 0
FullNghs: 1979.5 ave 2040 max 1895 min
Histogram: 1 0 0 0 1 0 0 0 0 2
Total # of neighbors = 7918
Ave neighs/atom = 15.4951
Neighbor list builds = 19
Dangerous builds = 0
Total wall time: 0:00:00

View File

@ -0,0 +1,92 @@
LAMMPS (4 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.SiC
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 edip/multi
pair_coeff * * SiC.edip Si C
Reading potential file SiC.edip with DATE: 2017-05-16
mass 1 28.085
mass 2 12.001
neighbor 1.0 bin
neigh_modify delay 1
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 1 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 3.94159
ghost atom cutoff = 3.94159
binsize = 1.97079, bins = 7 7 7
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair edip/multi, perpetual
attributes: full, newton on
pair build: full/bin/atomonly
stencil: full/bin/3d
bin: standard
Per MPI rank memory allocation (min/avg/max) = 2.692 | 2.692 | 2.692 Mbytes
Step Temp E_pair E_mol TotEng Press
0 0 -563.61621 0 -563.61621 -726147.34
10 4224.3601 -633.24829 0 -563.90103 -312355.55
20 4528.5661 -638.15183 0 -563.81071 -20091.291
30 4817.3654 -642.92111 0 -563.83905 106625.5
40 4619.4324 -639.6884 0 -563.85562 107180.42
50 4783.0025 -642.26961 0 -563.75166 75134.335
60 4525.145 -638.06177 0 -563.77681 71591.713
70 4685.2578 -640.72377 0 -563.8104 63956.042
80 4621.8393 -639.75912 0 -563.88682 18177.383
90 4834.7702 -643.34582 0 -563.97805 15282.823
100 4424.0589 -636.60208 0 -563.97656 47963.501
Loop time of 0.0552888 on 1 procs for 100 steps with 128 atoms
Performance: 156.270 ns/day, 0.154 hours/ns, 1808.685 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.051872 | 0.051872 | 0.051872 | 0.0 | 93.82
Neigh | 0.0023525 | 0.0023525 | 0.0023525 | 0.0 | 4.25
Comm | 0.0004518 | 0.0004518 | 0.0004518 | 0.0 | 0.82
Output | 0.00014806 | 0.00014806 | 0.00014806 | 0.0 | 0.27
Modify | 0.00024796 | 0.00024796 | 0.00024796 | 0.0 | 0.45
Other | | 0.0002165 | | | 0.39
Nlocal: 128 ave 128 max 128 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 473 ave 473 max 473 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 0 ave 0 max 0 min
Histogram: 1 0 0 0 0 0 0 0 0 0
FullNghs: 2376 ave 2376 max 2376 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 2376
Ave neighs/atom = 18.5625
Neighbor list builds = 11
Dangerous builds = 0
Total wall time: 0:00:00

View File

@ -0,0 +1,92 @@
LAMMPS (4 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.SiC
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 edip/multi
pair_coeff * * SiC.edip Si C
Reading potential file SiC.edip with DATE: 2017-05-16
mass 1 28.085
mass 2 12.001
neighbor 1.0 bin
neigh_modify delay 1
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 1 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 3.94159
ghost atom cutoff = 3.94159
binsize = 1.97079, bins = 7 7 7
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair edip/multi, perpetual
attributes: full, newton on
pair build: full/bin/atomonly
stencil: full/bin/3d
bin: standard
Per MPI rank memory allocation (min/avg/max) = 2.686 | 2.686 | 2.686 Mbytes
Step Temp E_pair E_mol TotEng Press
0 0 -563.61621 0 -563.61621 -726147.34
10 4224.3601 -633.24829 0 -563.90103 -312355.55
20 4528.5661 -638.15183 0 -563.81071 -20091.291
30 4817.3654 -642.92111 0 -563.83905 106625.5
40 4619.4324 -639.6884 0 -563.85562 107180.42
50 4783.0025 -642.26961 0 -563.75166 75134.335
60 4525.145 -638.06177 0 -563.77681 71591.713
70 4685.2578 -640.72377 0 -563.8104 63956.042
80 4621.8393 -639.75912 0 -563.88682 18177.383
90 4834.7702 -643.34582 0 -563.97805 15282.823
100 4424.0589 -636.60208 0 -563.97656 47963.501
Loop time of 0.020755 on 4 procs for 100 steps with 128 atoms
Performance: 416.285 ns/day, 0.058 hours/ns, 4818.118 timesteps/s
99.2% 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.011816 | 0.013825 | 0.016871 | 1.6 | 66.61
Neigh | 0.00061321 | 0.00066817 | 0.00074816 | 0.0 | 3.22
Comm | 0.0023363 | 0.0054012 | 0.0075014 | 2.7 | 26.02
Output | 0.00020909 | 0.00022268 | 0.00025558 | 0.0 | 1.07
Modify | 8.3208e-05 | 9.346e-05 | 0.00010395 | 0.0 | 0.45
Other | | 0.0005446 | | | 2.62
Nlocal: 32 ave 36 max 25 min
Histogram: 1 0 0 0 0 0 0 1 1 1
Nghost: 262.75 ave 273 max 255 min
Histogram: 2 0 0 0 0 0 0 1 0 1
Neighs: 0 ave 0 max 0 min
Histogram: 4 0 0 0 0 0 0 0 0 0
FullNghs: 594 ave 687 max 453 min
Histogram: 1 0 0 0 0 0 1 1 0 1
Total # of neighbors = 2376
Ave neighs/atom = 18.5625
Neighbor list builds = 11
Dangerous builds = 0
Total wall time: 0:00:00

2
src/.gitignore vendored
View File

@ -637,6 +637,8 @@
/pair_eam_fs_opt.h
/pair_edip.cpp
/pair_edip.h
/pair_edip_multi.cpp
/pair_edip_multi.h
/pair_eff_cut.cpp
/pair_eff_cut.h
/pair_eff_inline.h

View File

@ -798,6 +798,9 @@ void PairEDIP::coeff(int narg, char **arg)
}
}
if (nelements != 1)
error->all(FLERR,"Pair style edip only supports single element potentials");
// read potential file and initialize potential parameters
read_file(arg[2]);
@ -836,7 +839,7 @@ void PairEDIP::coeff(int narg, char **arg)
void PairEDIP::init_style()
{
if (force->newton_pair == 0)
error->all(FLERR,"Pair style EDIP requires newton pair on");
error->all(FLERR,"Pair style edip requires newton pair on");
// need a full neighbor list

View File

@ -0,0 +1,784 @@
/* ----------------------------------------------------------------------
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.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Environment Dependent Interatomic Potential
Contributing author: Chao Jiang
------------------------------------------------------------------------- */
#include <math.h>
#include <float.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "pair_edip_multi.h"
#include "atom.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "force.h"
#include "comm.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
#define MAXLINE 1024
#define DELTA 4
/* ---------------------------------------------------------------------- */
PairEDIPMulti::PairEDIPMulti(LAMMPS *lmp) : Pair(lmp)
{
single_enable = 0;
restartinfo = 0;
one_coeff = 1;
manybody_flag = 1;
nelements = 0;
elements = NULL;
nparams = maxparam = 0;
params = NULL;
elem2param = NULL;
}
/* ----------------------------------------------------------------------
check if allocated, since class can be destructed when incomplete
------------------------------------------------------------------------- */
PairEDIPMulti::~PairEDIPMulti()
{
if (elements)
for (int i = 0; i < nelements; i++) delete [] elements[i];
delete [] elements;
memory->destroy(params);
memory->destroy(elem2param);
if (allocated) {
memory->destroy(setflag);
memory->destroy(cutsq);
delete [] map;
//XXX deallocateGrids();
deallocatePreLoops();
}
}
/* ---------------------------------------------------------------------- */
void PairEDIPMulti::compute(int eflag, int vflag)
{
int i,j,k,ii,jj,kk,inum,jnum;
int itype,jtype,ktype,ijparam,ikparam,ijkparam;
double xtmp,ytmp,ztmp,evdwl;
int *ilist,*jlist,*numneigh,**firstneigh;
register int preForceCoord_counter;
double zeta_i;
double dzetair;
double fpair;
double costheta;
double dpairZ,dtripleZ;
// eflag != 0 means compute energy contributions in this step
// vflag != 0 means compute virial contributions in this step
evdwl = 0.0;
if (eflag || vflag) ev_setup(eflag,vflag);
else evflag = vflag_fdotr = 0;
double **x = atom->x;
double **f = atom->f;
int *type = atom->type;
int nlocal = atom->nlocal;
int newton_pair = force->newton_pair;
inum = list->inum;//total number of atoms in the cell
ilist = list->ilist;//list of atoms
numneigh = list->numneigh;//number of near neighbors
firstneigh = list->firstneigh;//list of neighbors
// loop over full neighbor list of my atoms
for (ii = 0; ii < inum; ii++) {
zeta_i = 0.0;
int numForceCoordPairs = 0;
i = ilist[ii];
itype = map[type[i]];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
// all the neighbors of atom i
jlist = firstneigh[i];
jnum = numneigh[i];
// pre-loop to compute environment coordination f(Z)
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
double delx, dely, delz, r_ij;
delx = x[j][0] - xtmp;
dely = x[j][1] - ytmp;
delz = x[j][2] - ztmp;
r_ij = delx * delx + dely * dely + delz * delz;
jtype = map[type[j]];
ijparam = elem2param[itype][jtype][jtype];
if (r_ij > params[ijparam].cutsq) continue;
r_ij = sqrt(r_ij);
// zeta and its derivative dZ/dr
if (r_ij < params[ijparam].cutoffC) zeta_i += 1.0;
else {
double f, fdr;
edip_fc(r_ij, &params[ijparam], f, fdr);
zeta_i += f;
dzetair = -fdr / r_ij;
preForceCoord_counter=numForceCoordPairs*5;
preForceCoord[preForceCoord_counter+0]=dzetair;
preForceCoord[preForceCoord_counter+1]=delx;
preForceCoord[preForceCoord_counter+2]=dely;
preForceCoord[preForceCoord_counter+3]=delz;
preForceCoord[preForceCoord_counter+4]=j;
numForceCoordPairs++;
}
}
// two-body interactions
dpairZ=0;
dtripleZ=0;
for (jj = 0; jj < jnum; jj++) {
double dr_ij[3], r_ij, f_ij[3];
j = jlist[jj];
j &= NEIGHMASK;
dr_ij[0] = x[j][0] - xtmp;
dr_ij[1] = x[j][1] - ytmp;
dr_ij[2] = x[j][2] - ztmp;
r_ij = dr_ij[0]*dr_ij[0] + dr_ij[1]*dr_ij[1] + dr_ij[2]*dr_ij[2];
jtype = map[type[j]];
ijparam = elem2param[itype][jtype][jtype];
if (r_ij > params[ijparam].cutsq) continue;
r_ij = sqrt(r_ij);
// potential energy and force
// since pair i-j is different from pair j-i, double counting is
// already considered in constructing the potential
double fdr, fdZ;
edip_pair(r_ij, zeta_i, &params[ijparam], evdwl, fdr, fdZ);
fpair = -fdr / r_ij;
dpairZ += fdZ;
f[i][0] -= fpair * dr_ij[0];
f[i][1] -= fpair * dr_ij[1];
f[i][2] -= fpair * dr_ij[2];
f[j][0] += fpair * dr_ij[0];
f[j][1] += fpair * dr_ij[1];
f[j][2] += fpair * dr_ij[2];
if (evflag) ev_tally(i, j, nlocal, newton_pair, evdwl, 0.0, fpair, -dr_ij[0], -dr_ij[1], -dr_ij[2]);
// three-body Forces
for (kk = jj + 1; kk < jnum; kk++) {
double dr_ik[3], r_ik, f_ik[3];
k = jlist[kk];
k &= NEIGHMASK;
ktype = map[type[k]];
ikparam = elem2param[itype][ktype][ktype];
ijkparam = elem2param[itype][jtype][ktype];
dr_ik[0] = x[k][0] - xtmp;
dr_ik[1] = x[k][1] - ytmp;
dr_ik[2] = x[k][2] - ztmp;
r_ik = dr_ik[0]*dr_ik[0] + dr_ik[1]*dr_ik[1] + dr_ik[2]*dr_ik[2];
if (r_ik > params[ikparam].cutsq) continue;
r_ik = sqrt(r_ik);
costheta=vec3_dot(dr_ij, dr_ik) / r_ij / r_ik;
double v1, v2, v3, v4, v5, v6, v7;
edip_fcut3(r_ij, &params[ijparam], v1, v2);
edip_fcut3(r_ik, &params[ikparam], v3, v4);
edip_h(costheta, zeta_i, &params[ijkparam], v5, v6, v7);
// potential energy and forces
evdwl = v1 * v3 * v5;
dtripleZ += v1 * v3 * v7;
double dri[3], drj[3], drk[3];
double dhl, dfr;
dhl = v1 * v3 * v6;
costheta_d(dr_ij, r_ij, dr_ik, r_ik, dri, drj, drk);
f_ij[0] = -dhl * drj[0];
f_ij[1] = -dhl * drj[1];
f_ij[2] = -dhl * drj[2];
f_ik[0] = -dhl * drk[0];
f_ik[1] = -dhl * drk[1];
f_ik[2] = -dhl * drk[2];
dfr = v2 * v3 * v5;
fpair = -dfr / r_ij;
f_ij[0] += fpair * dr_ij[0];
f_ij[1] += fpair * dr_ij[1];
f_ij[2] += fpair * dr_ij[2];
dfr = v1 * v4 * v5;
fpair = -dfr / r_ik;
f_ik[0] += fpair * dr_ik[0];
f_ik[1] += fpair * dr_ik[1];
f_ik[2] += fpair * dr_ik[2];
f[j][0] += f_ij[0];
f[j][1] += f_ij[1];
f[j][2] += f_ij[2];
f[k][0] += f_ik[0];
f[k][1] += f_ik[1];
f[k][2] += f_ik[2];
f[i][0] -= f_ij[0] + f_ik[0];
f[i][1] -= f_ij[1] + f_ik[1];
f[i][2] -= f_ij[2] + f_ik[2];
if (evflag) ev_tally3(i,j,k,evdwl,0.0,f_ij,f_ik,dr_ij,dr_ik);
}
}
// forces due to environment coordination f(Z)
for (int idx = 0; idx < numForceCoordPairs; idx++) {
double delx, dely, delz;
preForceCoord_counter = idx * 5;
dzetair = preForceCoord[preForceCoord_counter+0];
delx = preForceCoord[preForceCoord_counter+1];
dely = preForceCoord[preForceCoord_counter+2];
delz = preForceCoord[preForceCoord_counter+3];
j = static_cast<int> (preForceCoord[preForceCoord_counter+4]);
dzetair *= (dpairZ + dtripleZ);
f[j][0] += dzetair * delx;
f[j][1] += dzetair * dely;
f[j][2] += dzetair * delz;
f[i][0] -= dzetair * delx;
f[i][1] -= dzetair * dely;
f[i][2] -= dzetair * delz;
evdwl = 0.0;
if (evflag) ev_tally(i, j, nlocal, newton_pair, evdwl, 0.0, dzetair, -delx, -dely, -delz);
}
}
if (vflag_fdotr) virial_fdotr_compute();
}
double sqr(double x)
{
return x * x;
}
//pair Vij, partial derivatives dVij(r,Z)/dr and dVij(r,Z)/dZ
void PairEDIPMulti::edip_pair(double r, double z, Param *param, double &eng,
double &fdr, double &fZ)
{
double A = param->A;
double B = param->B;
double rho = param->rho;
double beta = param->beta;
double v1,v2,v3,v4;
v1 = pow(B / r, rho);
v2 = exp(-beta * z * z);
edip_fcut2(r, param, v3, v4);
eng = A * (v1 - v2) * v3;
fdr = A * (v1 - v2) * v4 + A * (-rho * v1 / r) * v3;
fZ = A * (2 * beta * z * v2) * v3;
}
//function fc(r) in calculating coordination Z and derivative fc'(r)
void PairEDIPMulti::edip_fc(double r, Param *param, double &f, double &fdr)
{
double a = param->cutoffA;
double c = param->cutoffC;
double alpha = param->alpha;
double x;
double v1, v2, v3;
if(r < c + 1E-6)
{
f=1.0;
fdr=0.0;
return;
}
if(r > a - 1E-6)
{
f=0.0;
fdr=0.0;
return;
}
x = (a - c) / (r - c);
v1 = x * x * x;
v2 = 1.0 / (1.0 - v1);
f = exp(alpha * v2);
fdr = (3.0 * x * v1 / (a - c)) * (-alpha * v2 * v2) * f;
}
//cut-off function for Vij and its derivative fcut2'(r)
void PairEDIPMulti::edip_fcut2(double r, Param *param, double &f, double &fdr)
{
double sigma = param->sigma;
double a = param->cutoffA;
double v1;
if(r > a - 1E-6)
{
f=0.0;
fdr=0.0;
return;
}
v1 = 1.0 / (r - a);
f = exp(sigma * v1);
fdr = (-sigma * v1 * v1) * f;
}
//function tau(Z) and its derivative tau'(Z)
void PairEDIPMulti::edip_tau(double z, Param *param, double &f, double &fdZ)
{
double u1 = param->u1;
double u2 = param->u2;
double u3 = param->u3;
double u4 = param->u4;
double v1, v2;
v1 = exp(-u4 * z);
v2 = exp(-2.0 * u4 * z);
f = u1 + u2 * u3 * v1 - u2 * v2;
fdZ = -u2 * u3 * u4 * v1 + 2.0 * u2 * u4 * v2;
}
//function h(l,Z) and its partial derivatives dh(l,Z)/dl and dh(l,Z)/dZ
void PairEDIPMulti::edip_h(double l, double z, Param *param, double &f,
double &fdl, double &fdZ)
{
double lambda = param->lambda;
double eta = param->eta;
double Q0 = param->Q0;
double mu = param->mu;
double Q, QdZ, Tau, TaudZ;
double u2, du2l, du2Z;
double v1, v2, v3;
//function Q(Z)
Q = Q0 * exp(-mu * z);
//derivative Q'(Z)
QdZ= -mu * Q;
edip_tau(z, param, Tau, TaudZ);
v1 = sqr(l + Tau);
u2 = Q * v1;
v2 = exp(-u2);
f = lambda * (1 - v2 + eta * u2);
//df/du2
v3 = lambda * (v2 + eta);
//du2/dl
du2l = Q * 2 * (l + Tau);
fdl = v3 * du2l;
//du2/dZ
du2Z = QdZ * v1 + Q * 2 * (l + Tau) * TaudZ;
fdZ = v3 * du2Z;
}
//cut-off function for Vijk and its derivative fcut3'(r)
void PairEDIPMulti::edip_fcut3(double r, Param *param, double &f, double &fdr)
{
double gamma = param->gamma;
double a = param->cutoffA;
double v1;
if(r > a - 1E-6)
{
f=0.0;
fdr=0.0;
return;
}
v1 = 1.0 / (r - a);
f = exp(gamma * v1);
fdr = (-gamma * v1 * v1) * f;
}
/* ----------------------------------------------------------------------
pre-calculated structures
------------------------------------------------------------------------- */
void PairEDIPMulti::allocatePreLoops(void)
{
int nthreads = comm->nthreads;
memory->create(preForceCoord,5*nthreads*leadDimInteractionList,"edip:preForceCoord");
}
/* ----------------------------------------------------------------------
deallocate preLoops
------------------------------------------------------------------------- */
void PairEDIPMulti::deallocatePreLoops(void)
{
memory->destroy(preForceCoord);
}
/* ---------------------------------------------------------------------- */
void PairEDIPMulti::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 PairEDIPMulti::settings(int narg, char **arg)
{
if (narg != 0) error->all(FLERR,"Illegal pair_style command");
}
/* ----------------------------------------------------------------------
set coeffs for one or more type pairs
------------------------------------------------------------------------- */
void PairEDIPMulti::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");
// insure I,J args are * *
if (strcmp(arg[0],"*") != 0 || strcmp(arg[1],"*") != 0)
error->all(FLERR,"Incorrect args for pair coefficients");
// 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++;
}
}
// read potential file and initialize potential parameters
read_file(arg[2]);
setup();
// 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
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;
count++;
}
if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
// allocate tables and internal structures
allocatePreLoops();
}
/* ----------------------------------------------------------------------
init specific to this pair style
------------------------------------------------------------------------- */
void PairEDIPMulti::init_style()
{
if (atom->tag_enable == 0)
error->all(FLERR,"Pair style edip/multi requires atom IDs");
if (force->newton_pair == 0)
error->all(FLERR,"Pair style edip/multi requires newton pair on");
// need a full neighbor list
int irequest = neighbor->request(this,instance_me);
neighbor->requests[irequest]->half = 0;
neighbor->requests[irequest]->full = 1;
}
/* ----------------------------------------------------------------------
init for one type pair i,j and corresponding j,i
------------------------------------------------------------------------- */
double PairEDIPMulti::init_one(int i, int j)
{
if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set");
return cutmax;
}
/* ---------------------------------------------------------------------- */
void PairEDIPMulti::read_file(char *file)
{
int params_per_line = 20;
char **words = new char*[params_per_line+1];
memory->sfree(params);
params = NULL;
nparams = maxparam = 0;
// open file on proc 0
FILE *fp;
if (comm->me == 0) {
fp = force->open_potential(file);
if (fp == NULL) {
char str[128];
sprintf(str,"Cannot open EDIP potential file %s",file);
error->one(FLERR,str);
}
}
// read each set of params from potential file
// one set of params can span multiple lines
// store params if all 3 element tags are in element list
int n,nwords,ielement,jelement,kelement;
char line[MAXLINE],*ptr;
int 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';
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 EDIP potential file");
// words = ptrs to all words in line
nwords = 0;
words[nwords++] = strtok(line," \t\n\r\f");
while ((words[nwords++] = strtok(NULL," \t\n\r\f"))) continue;
// ielement,jelement,kelement = 1st args
// if all 3 args are in element list, then parse this line
// else skip to next entry in file
for (ielement = 0; ielement < nelements; ielement++)
if (strcmp(words[0],elements[ielement]) == 0) break;
if (ielement == nelements) continue;
for (jelement = 0; jelement < nelements; jelement++)
if (strcmp(words[1],elements[jelement]) == 0) break;
if (jelement == nelements) continue;
for (kelement = 0; kelement < nelements; kelement++)
if (strcmp(words[2],elements[kelement]) == 0) break;
if (kelement == nelements) continue;
// load up parameter settings and error check their values
if (nparams == maxparam) {
maxparam += DELTA;
params = (Param *) memory->srealloc(params,maxparam*sizeof(Param),
"pair:params");
}
params[nparams].ielement = ielement;
params[nparams].jelement = jelement;
params[nparams].kelement = kelement;
params[nparams].A = atof(words[3]);
params[nparams].B = atof(words[4]);
params[nparams].cutoffA = atof(words[5]);
params[nparams].cutoffC = atof(words[6]);
params[nparams].alpha = atof(words[7]);
params[nparams].beta = atof(words[8]);
params[nparams].eta = atof(words[9]);
params[nparams].gamma = atof(words[10]);
params[nparams].lambda = atof(words[11]);
params[nparams].mu = atof(words[12]);
params[nparams].rho = atof(words[13]);
params[nparams].sigma = atof(words[14]);
params[nparams].Q0 = atof(words[15]);
params[nparams].u1 = atof(words[16]);
params[nparams].u2 = atof(words[17]);
params[nparams].u3 = atof(words[18]);
params[nparams].u4 = atof(words[19]);
if (params[nparams].A < 0.0 || params[nparams].B < 0.0 ||
params[nparams].cutoffA < 0.0 || params[nparams].cutoffC < 0.0 ||
params[nparams].alpha < 0.0 || params[nparams].beta < 0.0 ||
params[nparams].eta < 0.0 || params[nparams].gamma < 0.0 ||
params[nparams].lambda < 0.0 || params[nparams].mu < 0.0 ||
params[nparams].rho < 0.0 || params[nparams].sigma < 0.0)
error->all(FLERR,"Illegal EDIP parameter");
nparams++;
}
delete [] words;
}
/* ---------------------------------------------------------------------- */
void PairEDIPMulti::setup()
{
int i,j,k,m,n;
double rtmp;
// set elem2param for all triplet combinations
// must be a single exact match to lines read from file
// do not allow for ACB in place of ABC
memory->destroy(elem2param);
memory->create(elem2param,nelements,nelements,nelements,"pair:elem2param");
for (i = 0; i < nelements; i++)
for (j = 0; j < nelements; j++)
for (k = 0; k < nelements; k++) {
n = -1;
for (m = 0; m < nparams; m++) {
if (i == params[m].ielement && j == params[m].jelement &&
k == params[m].kelement) {
if (n >= 0) error->all(FLERR,"Potential file has duplicate entry");
n = m;
}
}
if (n < 0) error->all(FLERR,"Potential file is missing an entry");
elem2param[i][j][k] = n;
}
// set cutoff square
for (m = 0; m < nparams; m++) {
params[m].cutsq = params[m].cutoffA*params[m].cutoffA;
}
// set cutmax to max of all params
cutmax = 0.0;
for (m = 0; m < nparams; m++) {
rtmp = sqrt(params[m].cutsq);
if (rtmp > cutmax) cutmax = rtmp;
}
}

View File

@ -0,0 +1,113 @@
/* ----------------------------------------------------------------------
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(edip/multi,PairEDIPMulti)
#else
#ifndef LMP_PAIR_EDIP_MULTI_H
#define LMP_PAIR_EDIP_MULTI_H
#include "pair.h"
namespace LAMMPS_NS {
class PairEDIPMulti : public Pair {
public:
PairEDIPMulti(class LAMMPS *);
virtual ~PairEDIPMulti();
virtual void compute(int, int);
void settings(int, char **);
void coeff(int, char **);
double init_one(int, int);
void init_style();
protected:
struct Param {
double A, B;//coefficients for pair interaction I-J
double cutoffA;//cut-off distance for pair interaction I-J
double cutoffC;//lower cut-off distance for calculating Z_I
double alpha;//coefficient for calculating Z_I
double beta;//attractive term for pair I-J
double sigma;//cut-off coefficient for pair I-J
double rho;//pair I-J
double gamma;//coefficient for three-body interaction I-J-K
double eta, lambda;//coefficients for function h(l,Z)
double mu, Q0;//coefficients for function Q(Z)
double u1, u2, u3, u4;//coefficients for function tau(Z)
double cutsq;
int ielement,jelement,kelement;
};
double *preForceCoord;
double cutmax; // max cutoff for all elements
int nelements; // # of unique elements
char **elements; // names of unique elements
int ***elem2param; // mapping from element triplets to parameters
int *map; // mapping from atom types to elements
int nparams; // # of stored parameter sets
int maxparam; // max # of parameter sets
Param *params; // parameter set for an I-J-K interaction
// max number of interaction per atom for f(Z) environment potential
static const int leadDimInteractionList = 64;
void allocate();
void allocatePreLoops(void);
void deallocatePreLoops(void);
void read_file(char *);
void setup();
void edip_pair(double, double, Param *, double &, double &, double &);
void edip_fc(double, Param *, double &, double &);
void edip_fcut2(double, Param *, double &, double &);
void edip_tau(double, Param *, double &, double &);
void edip_h(double, double, Param *, double &, double &, double &);
void edip_fcut3(double, Param *, double &, double &);
double vec3_dot(double x[3], double y[3])
{
return x[0]*y[0] + x[1]*y[1] + x[2]*y[2];
}
void vec3_add(double k1, double x[3], double k2, double y[3], double *z)
{
z[0] = k1 * x[0] + k2 * y[0];
z[1] = k1 * x[1] + k2 * y[1];
z[2] = k1 * x[2] + k2 * y[2];
}
//dr_ij=r_j - r_i
//dr_ik=r_k - r_i
void costheta_d(double *dr_ij, double r_ij, double *dr_ik, double r_ik,
double *dri, double *drj, double *drk)
{
double costheta;
costheta = vec3_dot(dr_ij, dr_ik) / r_ij / r_ik;
vec3_add(1 / r_ij / r_ik, dr_ik, -costheta / r_ij / r_ij, dr_ij, drj);
vec3_add(1 / r_ij / r_ik, dr_ij, -costheta / r_ik / r_ik, dr_ik, drk);
vec3_add(-1, drj, -1, drk, dri);
}
};
}
#endif
#endif