diff --git a/examples/meam/dump.xyz b/examples/meam/dump.xyz new file mode 100644 index 0000000000..60e60ac7e7 --- /dev/null +++ b/examples/meam/dump.xyz @@ -0,0 +1,130 @@ +128 + Atoms. Timestep: 0 +2 -2.93785 -4.45926 -4.81092 +2 5.62221 -2.7335 -1.71576 +2 -2.66146 -5.54311 1.63537 +2 -5.43268 -4.61746 5.94523 +2 5.86792 -0.112054 -3.58394 +2 -3.71746 -0.662331 -0.371479 +2 -5.07247 -2.56716 4.41035 +2 -3.39514 0.934113 4.93107 +2 -5.43476 1.95238 -5.61809 +2 -4.58847 2.29045 -1.05977 +2 -5.90587 0.621241 2.01276 +2 -4.76807 0.196574 4.32678 +2 -5.42289 5.25697 -4.51629 +2 -5.2684 -5.91937 -2.86487 +2 -2.86109 1.04847 2.02991 +2 -4.07111 5.3133 3.80095 +2 -0.194715 -4.16777 -5.69509 +2 -2.98927 -3.16474 -1.61739 +2 -0.912931 -4.38191 -0.160186 +2 -2.45137 -5.24665 4.88829 +2 -2.888 -0.163345 -3.34011 +1 -4.67388 -1.38073 -2.29468 +2 -0.697395 -1.48853 0.600516 +1 -2.73922 -2.47748 0.238719 +2 -2.65513 -2.723 2.63503 +1 -3.46443 -4.60281 3.38178 +2 0.722761 -2.07094 2.92147 +1 -2.10006 -3.21313 5.72734 +2 -3.10576 2.32048 -2.27256 +1 -2.22988 0.716839 -1.31072 +2 -1.86983 1.40068 0.726511 +1 -4.11034 -0.709334 1.93418 +2 -0.350558 3.27072 -0.288066 +1 -3.40454 -1.4384 4.39035 +2 -3.09405 1.41325 -5.36355 +1 -4.45607 1.20729 -3.73102 +2 -2.6061 4.63735 -4.69039 +1 -3.34774 4.67681 -2.62847 +2 0.81217 4.86024 -4.67109 +1 -2.57569 3.37407 -0.213635 +2 -0.386798 5.87456 -2.11199 +1 -1.67662 1.33743 3.87415 +2 -0.877061 3.37359 4.3847 +1 -1.86093 3.31582 -5.97866 +1 -5.27323 -4.60733 -0.958175 +1 -2.78887 -5.69102 -0.792202 +1 -2.47172 4.58019 2.50832 +1 -3.88199 5.84566 -5.75634 +2 2.23148 -2.77292 -5.23569 +2 0.298198 -3.13853 -3.16082 +2 2.88108 -3.46587 -0.58232 +2 0.250962 -5.75952 2.73898 +2 -0.293412 -0.802943 -3.36985 +1 -1.00757 -2.04819 -1.94193 +2 2.07291 1.49224 -2.38981 +1 1.11109 -3.20042 0.949108 +2 1.67743 -0.790186 2.51588 +1 -0.83423 -4.33425 2.09715 +2 3.27474 -1.31079 4.78847 +1 1.71262 -3.36915 4.5581 +2 0.47706 1.7769 -5.33399 +1 0.294439 0.589278 -2.20301 +2 2.20393 3.15576 -2.02768 +1 -0.0404494 0.476782 1.03964 +2 1.13959 2.37634 2.3481 +1 -0.973837 -1.63252 3.75386 +2 -0.3292 0.299699 5.27708 +1 -1.61856 -0.396427 -5.17712 +2 2.59999 -5.19777 5.82307 +1 -1.62707 2.32109 -3.62999 +2 3.65327 4.92826 -5.43193 +1 0.0788934 4.0241 -2.50115 +2 2.85565 2.61687 2.11255 +1 0.973899 2.62554 4.34121 +2 3.74529 3.45214 4.59464 +1 2.08052 4.7039 5.32803 +1 -1.03242 -5.8155 -4.32658 +1 0.762244 -4.36316 -1.31566 +1 0.326368 3.99374 1.61723 +1 -0.435011 -5.79971 4.59591 +2 3.91611 -4.60528 -3.31917 +2 1.92407 5.73451 -1.97543 +2 -5.97945 -4.23694 1.86465 +2 4.334 -4.48452 5.37374 +2 2.27555 -0.632774 -5.79318 +1 1.87282 -1.55049 -3.456 +2 3.45581 -1.10541 -1.83331 +1 4.37882 -1.94665 -0.328464 +2 2.59992 -3.7549 2.57406 +1 3.99839 -4.48566 1.19687 +2 -5.72956 -2.14757 -5.99636 +1 4.26641 -2.6989 -5.80055 +2 4.52547 2.29068 -3.47658 +1 2.36031 1.34164 -4.41738 +2 4.77671 1.40612 -0.752462 +1 1.80727 -0.783597 -0.458199 +2 4.4745 0.373622 2.10683 +1 3.60812 -1.73157 2.40191 +2 4.62814 -0.286541 4.47565 +1 1.79752 0.289353 4.23308 +2 5.83415 4.49865 -5.96645 +1 3.24013 4.16552 -3.507 +2 4.87203 4.871 -2.33644 +1 3.55265 1.22628 0.692683 +2 -5.81733 4.54205 1.55789 +1 3.96832 1.54411 3.82844 +2 -5.53493 1.9067 3.75041 +1 4.47286 2.64156 -5.59528 +1 1.70009 -4.81154 -4.19539 +1 1.72215 4.18784 -0.371268 +1 3.92182 4.59356 1.32634 +1 3.13102 -5.89225 3.60012 +1 4.75587 -2.28778 -3.47421 +1 -5.50503 -2.70274 0.874887 +1 5.84186 -4.60644 3.87141 +1 -4.75169 -3.1692 -4.40998 +1 3.9405 0.71887 -2.28988 +1 -5.68697 0.204238 -0.191674 +1 5.89496 -1.24226 3.12013 +1 5.96758 -0.0712572 5.8964 +1 -5.62085 3.36 -2.94935 +1 5.20653 3.45179 -0.380089 +1 -4.69945 2.54896 1.82974 +1 -4.07584 3.07262 5.0648 +1 4.15876 -5.08968 -1.14435 +1 -4.69638 -5.74298 1.13578 +1 5.59942 4.6887 3.59483 +1 5.09884 -5.37744 -4.90513 diff --git a/examples/meam/in.meam b/examples/meam/in.meam index b4463be365..550d2fb4dc 100644 --- a/examples/meam/in.meam +++ b/examples/meam/in.meam @@ -14,9 +14,10 @@ neighbor 0.3 bin neigh_modify delay 10 fix 1 all nve -thermo 10 +thermo 1 timestep 0.001 +dump 1 all xyz 100 dump.xyz #dump 1 all atom 50 dump.meam #dump 2 all image 10 image.*.jpg element element & @@ -27,4 +28,4 @@ timestep 0.001 # axes yes 0.8 0.02 view 60 -30 #dump_modify 3 pad 3 element Si C -run 100 +run 10 diff --git a/examples/meam/msmeam/test.dump b/examples/meam/msmeam/test.dump index 476ebfeb75..45d17a5228 100644 --- a/examples/meam/msmeam/test.dump +++ b/examples/meam/msmeam/test.dump @@ -7,9 +7,9 @@ ITEM: BOX BOUNDS pp pp pp -4.0000000000000000e+00 4.0000000000000000e+00 -4.0000000000000000e+00 4.0000000000000000e+00 ITEM: ATOMS id x y z fx fy fz c_pot_energy c_stress[1] c_stress[2] c_stress[3] c_stress[4] c_stress[5] c_stress[6] -1 0 0 0 -132.308 -88.6953 0 22.6425 -2.1528e+08 -1.63302e+08 -0 -2.06499e+07 -0 -0 -2 2.2 0 0 121.138 -0.446613 0 14.8521 -2.12596e+08 -0 -0 386632 -0 -0 -3 0.3 2.3 0 11.1693 89.1419 0 8.67443 -2.68428e+06 -1.63302e+08 -0 -2.10365e+07 -0 -0 +1 0 0 0 -131.977 -88.3322 0 22.9236 -2.14786e+08 -1.62719e+08 -0 -2.05378e+07 -0 -0 +2 2.2 0 0 120.857 -0.482171 0 14.7745 -2.12113e+08 -0 -0 403352 -0 -0 +3 0.3 2.3 0 11.1201 88.8144 0 8.61773 -2.67245e+06 -1.62719e+08 -0 -2.09411e+07 -0 -0 ITEM: TIMESTEP 1 ITEM: NUMBER OF ATOMS @@ -19,6 +19,6 @@ ITEM: BOX BOUNDS pp pp pp -4.0000000000000000e+00 4.0000000000000000e+00 -4.0000000000000000e+00 4.0000000000000000e+00 ITEM: ATOMS id x y z fx fy fz c_pot_energy c_stress[1] c_stress[2] c_stress[3] c_stress[4] c_stress[5] c_stress[6] -1 0 0 0 -132.308 -88.6953 0 22.6425 -2.1528e+08 -1.63302e+08 -0 -2.06499e+07 -0 -0 -2 2.2 0 0 121.138 -0.446613 0 14.8521 -2.12596e+08 -0 -0 386632 -0 -0 -3 0.3 2.3 0 11.1693 89.1419 0 8.67443 -2.68428e+06 -1.63302e+08 -0 -2.10365e+07 -0 -0 +1 0 0 0 -131.977 -88.3322 0 22.9236 -2.14786e+08 -1.62719e+08 -0 -2.05378e+07 -0 -0 +2 2.2 0 0 120.857 -0.482171 0 14.7745 -2.12113e+08 -0 -0 403352 -0 -0 +3 0.3 2.3 0 11.1201 88.8144 0 8.61773 -2.67245e+06 -1.62719e+08 -0 -2.09411e+07 -0 -0 diff --git a/examples/meam/msmeam/test.log b/examples/meam/msmeam/test.log index 83e0d92f70..813acd55da 100644 --- a/examples/meam/msmeam/test.log +++ b/examples/meam/msmeam/test.log @@ -83,24 +83,24 @@ Neighbor list info ... pair build: halffull/newton stencil: none bin: none -Per MPI rank memory allocation (min/avg/max) = 10.69 | 10.69 | 10.69 Mbytes +Per MPI rank memory allocation (min/avg/max) = 10.71 | 10.71 | 10.71 Mbytes Step Temp TotEng Press Pxx Pyy Pzz Pxy Pxz Pyz Lx Ly Lz Volume c_eatoms - 0 0 15.389692 492946.21 840939 637899.63 0 80663.676 0 0 8 8 8 512 15.389692 - 1 0 15.389692 492946.21 840939 637899.63 0 80663.676 0 0 8 8 8 512 15.389692 -Loop time of 7.9076e-05 on 1 procs for 1 steps with 3 atoms + 0 0 15.438614 491542.52 839006.02 635621.55 0 80225.587 0 0 8 8 8 512 15.438614 + 1 0 15.438614 491542.52 839006.02 635621.55 0 80225.587 0 0 8 8 8 512 15.438614 +Loop time of 6.7095e-05 on 1 procs for 1 steps with 3 atoms -Performance: 1092.620 ns/day, 0.022 hours/ns, 12646.062 timesteps/s -67.0% CPU use with 1 MPI tasks x no OpenMP threads +Performance: 1287.726 ns/day, 0.019 hours/ns, 14904.240 timesteps/s +77.5% CPU use with 1 MPI tasks x no OpenMP threads MPI task timing breakdown: Section | min time | avg time | max time |%varavg| %total --------------------------------------------------------------- -Pair | 2.4076e-05 | 2.4076e-05 | 2.4076e-05 | 0.0 | 30.45 +Pair | 2.7954e-05 | 2.7954e-05 | 2.7954e-05 | 0.0 | 41.66 Neigh | 0 | 0 | 0 | 0.0 | 0.00 -Comm | 3.204e-06 | 3.204e-06 | 3.204e-06 | 0.0 | 4.05 -Output | 4.726e-05 | 4.726e-05 | 4.726e-05 | 0.0 | 59.77 -Modify | 5.97e-07 | 5.97e-07 | 5.97e-07 | 0.0 | 0.75 -Other | | 3.939e-06 | | | 4.98 +Comm | 3.164e-06 | 3.164e-06 | 3.164e-06 | 0.0 | 4.72 +Output | 3.1658e-05 | 3.1658e-05 | 3.1658e-05 | 0.0 | 47.18 +Modify | 1.712e-06 | 1.712e-06 | 1.712e-06 | 0.0 | 2.55 +Other | | 2.607e-06 | | | 3.89 Nlocal: 3 ave 3 max 3 min Histogram: 1 0 0 0 0 0 0 0 0 0 diff --git a/examples/meam/msmeam/test.out b/examples/meam/msmeam/test.out new file mode 100644 index 0000000000..72592b4673 --- /dev/null +++ b/examples/meam/msmeam/test.out @@ -0,0 +1,294 @@ +LAMMPS (23 Jun 2022) +log test.log +# Test of MEAM potential for HGa + +# ------------------------ INITIALIZATION ---------------------------- +units metal +dimension 3 +boundary p p p +atom_style atomic +variable latparam equal 4.646 +variable ncell equal 3 + +# ----------------------- ATOM DEFINITION ---------------------------- +region box block -4 4 -4 4 -4 4 +create_box 2 box +Created orthogonal box = (-4 -4 -4) to (4 4 4) + 1 by 1 by 1 MPI processor grid + +# + +include potential.mod +# NOTE: This script can be modified for different pair styles +# See in.elastic for more info. + +variable Pu string H +print "potential chosen ${Pu}" +potential chosen H +# Choose potential +pair_style meam ms # the `ms` flag turns on MS-MEAM +print "we just executed" +we just executed + +pair_coeff * * library.MSmeam ${Pu} Ga4 HGaMS.meam ${Pu} Ga4 +pair_coeff * * library.MSmeam H Ga4 HGaMS.meam ${Pu} Ga4 +pair_coeff * * library.MSmeam H Ga4 HGaMS.meam H Ga4 +Reading MEAM library file library.MSmeam with DATE: 2018-09-22 +Tokenizer::next_double() current = 1.0 +Tokenizer::next_int() current = 1 +Tokenizer::next_double() current = 1.0079 +Tokenizer::next_double() current = 2.960 +Tokenizer::next_double() current = 2.960 +Tokenizer::next_double() current = 3.0 +Tokenizer::next_double() current = 1.0 +Tokenizer::next_double() current = 1.0 +Tokenizer::next_double() current = 1.0 +Tokenizer::next_double() current = 3.0 +Tokenizer::next_double() current = 1.0 +Tokenizer::next_double() current = 0.741 +Tokenizer::next_double() current = 2.235 +Tokenizer::next_double() current = 2.50 +Tokenizer::next_double() current = 1.0 +Tokenizer::next_double() current = 0.44721 +Tokenizer::next_double() current = 0.0 +Tokenizer::next_double() current = 0.00 +Tokenizer::next_double() current = 0.0 +Tokenizer::next_double() current = 0.31623 +Tokenizer::next_double() current = 0 +Tokenizer::next_double() current = 6.70 +Tokenizer::next_int() current = 0 +Tokenizer::next_double() current = 12.0 +Tokenizer::next_int() current = 31 +Tokenizer::next_double() current = 69.723 +Tokenizer::next_double() current = 4.42 +Tokenizer::next_double() current = 4.80 +Tokenizer::next_double() current = 3.10 +Tokenizer::next_double() current = 6.00 +Tokenizer::next_double() current = 0.00 +Tokenizer::next_double() current = 0.0 +Tokenizer::next_double() current = 0.0 +Tokenizer::next_double() current = 0.5 +Tokenizer::next_double() current = 4.247 +Tokenizer::next_double() current = 2.897 +Tokenizer::next_double() current = 0.97 +Tokenizer::next_double() current = 1.0 +Tokenizer::next_double() current = 1.649 +Tokenizer::next_double() current = 1.435 +Tokenizer::next_double() current = 0.00 +Tokenizer::next_double() current = 0.0 +Tokenizer::next_double() current = 0.0 +Tokenizer::next_double() current = 2.0 +Tokenizer::next_double() current = 0.70 +Tokenizer::next_int() current = 0 +^^^^^ Filling tm +^^^^^ Filling tm +Tokenizer::next_double() current = 1 +Tokenizer::next_double() current = 1 +Tokenizer::next_double() current = 0 +Tokenizer::next_double() current = 1 +Tokenizer::next_double() current = 5.9 +Tokenizer::next_int() current = 1 +Tokenizer::next_int() current = 1 +Tokenizer::next_double() current = 0.460 +Tokenizer::next_int() current = 1 +Tokenizer::next_int() current = 1 +Tokenizer::next_double() current = 0.460 +Tokenizer::next_int() current = 1 +Tokenizer::next_int() current = 1 +Tokenizer::next_int() current = 1 +Tokenizer::next_double() current = 1.3 +Tokenizer::next_int() current = 1 +Tokenizer::next_int() current = 1 +Tokenizer::next_int() current = 1 +Tokenizer::next_double() current = 2.80 +Tokenizer::next_int() current = 1 +Tokenizer::next_int() current = 1 +Tokenizer::next_double() current = 1 +Tokenizer::next_int() current = 2 +Tokenizer::next_double() current = 0.6 +Tokenizer::next_int() current = 2 +Tokenizer::next_int() current = 2 +Tokenizer::next_double() current = 0.097 +Tokenizer::next_int() current = 2 +Tokenizer::next_int() current = 2 +Tokenizer::next_double() current = 0.097 +Tokenizer::next_int() current = 2 +Tokenizer::next_int() current = 2 +Tokenizer::next_double() current = 1 +Tokenizer::next_int() current = 1 +Tokenizer::next_int() current = 2 +Tokenizer::next_double() current = 0.300 +Tokenizer::next_int() current = 1 +Tokenizer::next_int() current = 2 +Tokenizer::next_double() current = 0.300 +Tokenizer::next_int() current = 1 +Tokenizer::next_int() current = 2 +Tokenizer::next_int() current = 1 +Tokenizer::next_int() current = 2 +Tokenizer::next_double() current = 3.19 +Tokenizer::next_int() current = 1 +Tokenizer::next_int() current = 2 +Tokenizer::next_double() current = -0.48 +Tokenizer::next_int() current = 1 +Tokenizer::next_int() current = 2 +Tokenizer::next_double() current = 6.6 +Tokenizer::next_int() current = 1 +Tokenizer::next_int() current = 1 +Tokenizer::next_int() current = 2 +Tokenizer::next_double() current = 2.0 +Tokenizer::next_int() current = 2 +Tokenizer::next_int() current = 1 +Tokenizer::next_int() current = 2 +Tokenizer::next_double() current = 2.0 +Tokenizer::next_int() current = 1 +Tokenizer::next_int() current = 2 +Tokenizer::next_int() current = 1 +Tokenizer::next_double() current = 2.0 +Tokenizer::next_int() current = 2 +Tokenizer::next_int() current = 2 +Tokenizer::next_int() current = 1 +Tokenizer::next_double() current = 1.4 +Tokenizer::next_int() current = 1 +Tokenizer::next_int() current = 2 +Tokenizer::next_int() current = 2 +Tokenizer::next_double() current = 1.4 +Tokenizer::next_int() current = 1 +Tokenizer::next_int() current = 1 +Tokenizer::next_int() current = 2 +Tokenizer::next_double() current = 1.4 +Tokenizer::next_int() current = 1 +Tokenizer::next_int() current = 2 +Tokenizer::next_double() current = 1 +# Setup neighbor style +neighbor 1.0 nsq +neigh_modify once no every 1 delay 0 check yes + +# Setup minimization style +variable dmax equal 1.0e-2 +min_style cg +min_modify dmax ${dmax} line quadratic +min_modify dmax 0.01 line quadratic +compute eng all pe/atom +compute eatoms all reduce sum c_eng + +# Setup output +thermo 100 +thermo_style custom step temp etotal press pxx pyy pzz pxy pxz pyz lx ly lz vol c_eatoms +thermo_modify norm yes +create_atoms 1 single 0 0 0 units box +Created 1 atoms + using box units in orthogonal box = (-4 -4 -4) to (4 4 4) + create_atoms CPU = 0.000 seconds +create_atoms 2 single 2.2 0 0 units box +Created 1 atoms + using box units in orthogonal box = (-4 -4 -4) to (4 4 4) + create_atoms CPU = 0.000 seconds +create_atoms 2 single 0.3 2.3 0 units box +Created 1 atoms + using box units in orthogonal box = (-4 -4 -4) to (4 4 4) + create_atoms CPU = 0.000 seconds +# ---------- Define Settings --------------------- +variable teng equal "c_eatoms" +compute pot_energy all pe/atom +compute stress all stress/atom NULL +dump 1 all custom 1 test.dump id x y z fx fy fz c_pot_energy c_stress[1] c_stress[2] c_stress[3] c_stress[4] c_stress[5] c_stress[6] +run 1 +WARNING: No fixes with time integration, atoms won't move (src/verlet.cpp:60) +Neighbor list info ... + update every 1 steps, delay 0 steps, check yes + max neighbors/atom: 2000, page size: 100000 + master list distance cutoff = 6.9 + ghost atom cutoff = 6.9 + 2 neighbor lists, perpetual/occasional/extra = 2 0 0 + (1) pair meam, perpetual + attributes: full, newton on + pair build: full/nsq + stencil: none + bin: none + (2) pair meam, perpetual, half/full from (1) + attributes: half, newton on + pair build: halffull/newton + stencil: none + bin: none +Setting up Verlet run ... + Unit style : metal + Current step : 0 + Time step : 0.001 +--- compute derivatives of energy wrt ... +phip sij -116.571511 1.000000 +drhodr -1.327736 -0.011564 +fp 4.976162 -13.618158 +dUdrij -123.021058 +--- compute derivatives of energy wrt ... +phip sij -86.757953 1.000000 +drhodr -0.821855 -0.007170 +fp 4.976162 -14.961215 +dUdrij -90.740359 +--- compute derivatives of energy wrt ... +phip sij 0.129370 1.000000 +drhodr -0.000000 -0.002189 +fp -13.618158 4.976162 +dUdrij 0.118475 +--- compute derivatives of energy wrt ... +phip sij 0.140206 1.000000 +drhodr -0.000000 -0.002956 +fp -14.961215 4.976162 +dUdrij 0.125499 +Per MPI rank memory allocation (min/avg/max) = 10.71 | 10.71 | 10.71 Mbytes + Step Temp TotEng Press Pxx Pyy Pzz Pxy Pxz Pyz Lx Ly Lz Volume c_eatoms + 0 0 15.438614 491542.52 839006.02 635621.55 0 80225.587 0 0 8 8 8 512 15.438614 +--- compute derivatives of energy wrt ... +phip sij -116.571511 1.000000 +drhodr -1.327736 -0.011564 +fp 4.976162 -13.618158 +dUdrij -123.021058 +--- compute derivatives of energy wrt ... +phip sij -86.757953 1.000000 +drhodr -0.821855 -0.007170 +fp 4.976162 -14.961215 +dUdrij -90.740359 +--- compute derivatives of energy wrt ... +phip sij 0.129370 1.000000 +drhodr -0.000000 -0.002189 +fp -13.618158 4.976162 +dUdrij 0.118475 +--- compute derivatives of energy wrt ... +phip sij 0.140206 1.000000 +drhodr -0.000000 -0.002956 +fp -14.961215 4.976162 +dUdrij 0.125499 + 1 0 15.438614 491542.52 839006.02 635621.55 0 80225.587 0 0 8 8 8 512 15.438614 +Loop time of 6.7095e-05 on 1 procs for 1 steps with 3 atoms + +Performance: 1287.726 ns/day, 0.019 hours/ns, 14904.240 timesteps/s +77.5% CPU use with 1 MPI tasks x no OpenMP threads + +MPI task timing breakdown: +Section | min time | avg time | max time |%varavg| %total +--------------------------------------------------------------- +Pair | 2.7954e-05 | 2.7954e-05 | 2.7954e-05 | 0.0 | 41.66 +Neigh | 0 | 0 | 0 | 0.0 | 0.00 +Comm | 3.164e-06 | 3.164e-06 | 3.164e-06 | 0.0 | 4.72 +Output | 3.1658e-05 | 3.1658e-05 | 3.1658e-05 | 0.0 | 47.18 +Modify | 1.712e-06 | 1.712e-06 | 1.712e-06 | 0.0 | 2.55 +Other | | 2.607e-06 | | | 3.89 + +Nlocal: 3 ave 3 max 3 min +Histogram: 1 0 0 0 0 0 0 0 0 0 +Nghost: 78 ave 78 max 78 min +Histogram: 1 0 0 0 0 0 0 0 0 0 +Neighs: 7 ave 7 max 7 min +Histogram: 1 0 0 0 0 0 0 0 0 0 +FullNghs: 14 ave 14 max 14 min +Histogram: 1 0 0 0 0 0 0 0 0 0 + +Total # of neighbors = 14 +Ave neighs/atom = 4.6666667 +Neighbor list builds = 0 +Dangerous builds = 0 +write_data test.data +System init for write_data ... + +print "All done!" +All done! +Total wall time: 0:00:00 diff --git a/src/MEAM/meam_dens_final.cpp b/src/MEAM/meam_dens_final.cpp index 234694473a..201b27cbeb 100644 --- a/src/MEAM/meam_dens_final.cpp +++ b/src/MEAM/meam_dens_final.cpp @@ -25,8 +25,8 @@ MEAM::meam_dens_final(int nlocal, int eflag_either, int eflag_global, int eflag_ double denom, rho_bkgd, Fl; double scaleii; - printf("meam_dens_final arho3m\n"); - printf("%f\n", arho3m[0][0]); + //printf("meam_dens_final arho3m\n"); + //printf("%f\n", arho3m[0][0]); // Complete the calculation of density @@ -52,14 +52,14 @@ MEAM::meam_dens_final(int nlocal, int eflag_either, int eflag_global, int eflag_ for (m = 0; m < 10; m++) { rho3[i] = rho3[i] + this->v3D[m] * (arho3[i][m] * arho3[i][m] - arho3m[i][m] * arho3m[i][m]); - printf("arho loop %f %f\n", arho3[i][m], arho3m[i][m]); + //printf("arho loop %f %f\n", arho3[i][m], arho3m[i][m]); } // Compared to Greg's original, all of the t weightings are already accounted for. // Code block for t_ave removed. - printf("gam rho %f %f %f\n", rho1[i], rho2[i], rho3[i]); + //printf("gam rho %f %f %f\n", rho1[i], rho2[i], rho3[i]); gamma[i] = rho1[i] + rho2[i] + rho3[i]; - printf("gam check 0 %f\n", gamma[i]); + //printf("gam check 0 %f\n", gamma[i]); /* if (rho0[i] > 0.0) { @@ -80,11 +80,11 @@ MEAM::meam_dens_final(int nlocal, int eflag_either, int eflag_global, int eflag_ gamma[i] = t_ave[i][0] * rho1[i] + t_ave[i][1] * rho2[i] + t_ave[i][2] * rho3[i]; */ - printf("rho0 %f\n", rho0[i]); + //printf("rho0 %f\n", rho0[i]); if (rho0[i] > 0.0) { gamma[i] = gamma[i] / (rho0[i] * rho0[i]); } - printf("gam check1 %f\n", gamma[i]); + //printf("gam check1 %f\n", gamma[i]); Z = get_Zij(this->lattce_meam[elti][elti]); @@ -97,7 +97,7 @@ MEAM::meam_dens_final(int nlocal, int eflag_either, int eflag_global, int eflag_ if (this->ibar_meam[elti] <= 0) { Gbar = 1.0; dGbar = 0.0; - printf("ibar meam <= 0\n"); + //printf("ibar meam <= 0\n"); } else { if (this->mix_ref_t == 1) { gam = (t_ave[i][0] * shp[0] + t_ave[i][1] * shp[1] + t_ave[i][2] * shp[2]) / (Z * Z); @@ -106,7 +106,7 @@ MEAM::meam_dens_final(int nlocal, int eflag_either, int eflag_global, int eflag_ (Z * Z); } Gbar = G_gam(gam, this->ibar_meam[elti], errorflag); - printf("ibar meam > 0\n"); + //printf("ibar meam > 0\n"); } rho[i] = rho0[i] * G; @@ -130,7 +130,7 @@ MEAM::meam_dens_final(int nlocal, int eflag_either, int eflag_global, int eflag_ rhob = rho[i] / rho_bkgd; denom = 1.0 / rho_bkgd; - printf("dgam rhob denom %f %f\n", rhob, denom); + //printf("dgam rhob denom %f %f\n", rhob, denom); G = dG_gam(gamma[i], this->ibar_meam[elti], dG); diff --git a/src/MEAM/meam_dens_init.cpp b/src/MEAM/meam_dens_init.cpp index 802f811071..61fc6e8bf0 100644 --- a/src/MEAM/meam_dens_init.cpp +++ b/src/MEAM/meam_dens_init.cpp @@ -90,7 +90,7 @@ MEAM::meam_dens_setup(int atom_nmax, int nall, int n_neigh) } // zero out local arrays - printf("^^^ zero out local arrays in mean_dens_init\n"); + //printf("^^^ zero out local arrays in mean_dens_init\n"); for (i = 0; i < nall; i++) { rho0[i] = 0.0; arho2b[i] = 0.0; @@ -119,7 +119,7 @@ MEAM::meam_dens_setup(int atom_nmax, int nall, int n_neigh) tsq_ave[i][0] = tsq_ave[i][1] = tsq_ave[i][2] = 0.0; } - printf("^^^ arho3m %f\n", arho3m[0][0]); + //printf("^^^ arho3m %f\n", arho3m[0][0]); } void @@ -361,9 +361,11 @@ MEAM::calc_rho1(int i, int /*ntype*/, int* type, int* fmap, double** x, int numn rhoa2i = rhoa2i * this->t2_meam[elti]; rhoa3i = rhoa3i * this->t3_meam[elti]; } + /* printf("rhos1 %d %d %f %f %f %f\n", i, j, rhoa1i, rhoa1j, rhoa1mi, rhoa1mj); printf("rhos2 %d %d %f %f %f %f\n", i, j, rhoa2i, rhoa2j, rhoa2mi, rhoa2mj); printf("rhos3 %d %d %f %f %f %f\n", i, j, rhoa3i, rhoa3j, rhoa3mi, rhoa3mj); + */ rho0[i] = rho0[i] + rhoa0j; rho0[j] = rho0[j] + rhoa0i; // For ialloy = 2, use single-element value (not average) @@ -405,7 +407,7 @@ MEAM::calc_rho1(int i, int /*ntype*/, int* type, int* fmap, double** x, int numn A2mi = rhoa2mi/rij2; A3mi = rhoa3mi/(rij2*rij); } - printf("arho2mb i j %d %d %f %f\n", i, j, arho2mb[i], arho2mb[j]); + //printf("arho2mb i j %d %d %f %f\n", i, j, arho2mb[i], arho2mb[j]); for (m = 0; m < 3; m++) { arho1[i][m] = arho1[i][m] + A1j * delij[m]; arho1[j][m] = arho1[j][m] - A1i * delij[m]; @@ -417,8 +419,8 @@ MEAM::calc_rho1(int i, int /*ntype*/, int* type, int* fmap, double** x, int numn arho3mb[i][m] = arho3mb[i][m] + rhoa3mj*delij[m] / rij; arho3mb[j][m] = arho3mb[j][m] - rhoa3mi*delij[m] / rij; } - printf("arho1m %f %f\n", arho1m[i][m], arho1m[j][m]); - printf("arho3mb %f %f\n", arho3mb[i][m], arho3mb[j][m]); + //printf("arho1m %f %f\n", arho1m[i][m], arho1m[j][m]); + //printf("arho3mb %f %f\n", arho3mb[i][m], arho3mb[j][m]); for (n = m; n < 3; n++) { arho2[i][nv2] = arho2[i][nv2] + A2j * delij[m] * delij[n]; arho2[j][nv2] = arho2[j][nv2] + A2i * delij[m] * delij[n]; @@ -427,7 +429,7 @@ MEAM::calc_rho1(int i, int /*ntype*/, int* type, int* fmap, double** x, int numn arho2m[j][nv2] = arho2m[j][nv2] + A2mi*delij[m] * delij[n]; } //printf("delij %f %f\n", delij[m], delij[n]); - printf("arho2m %d %d %f %f\n", nv2, m, arho2m[i][nv2], arho2m[j][nv2]); + //printf("arho2m %d %d %f %f\n", nv2, m, arho2m[i][nv2], arho2m[j][nv2]); nv2 = nv2 + 1; for (p = n; p < 3; p++) { arho3[i][nv3] = arho3[i][nv3] + A3j * delij[m] * delij[n] * delij[p]; @@ -437,7 +439,7 @@ MEAM::calc_rho1(int i, int /*ntype*/, int* type, int* fmap, double** x, int numn arho3m[j][nv3] = arho3m[j][nv3] - A3mi*delij[m]*delij[n]*delij[p]; } //arho3m[0][0]=500.0; - printf("arho3m %d %f %f\n", p, arho3m[i][nv3], arho3m[j][nv3]); + //printf("arho3m %d %f %f\n", p, arho3m[i][nv3], arho3m[j][nv3]); nv3 = nv3 + 1; } } diff --git a/src/MEAM/meam_force.cpp b/src/MEAM/meam_force.cpp index bca423f563..172f61a2dc 100644 --- a/src/MEAM/meam_force.cpp +++ b/src/MEAM/meam_force.cpp @@ -75,8 +75,8 @@ MEAM::meam_force(int i, int eflag_global, int eflag_atom, int vflag_global, int double drho3mdr1, drho3mdr2, drho3mds1, drho3mds2; double drho3mdrm1[3], drho3mdrm2[3]; - printf("----- begin meam_force -----\n"); - printf("forces f %d %f\n", i, f[0][0]); + //printf("----- begin meam_force -----\n"); + //printf("forces f %d %f\n", i, f[0][0]); third = 1.0 / 3.0; sixth = 1.0 / 6.0; @@ -89,7 +89,7 @@ MEAM::meam_force(int i, int eflag_global, int eflag_atom, int vflag_global, int xitmp = x[i][0]; yitmp = x[i][1]; zitmp = x[i][2]; - printf("xtmp ytmp ztmp: %f %f %f\n", xitmp, yitmp, zitmp); + //printf("xtmp ytmp ztmp: %f %f %f\n", xitmp, yitmp, zitmp); // Treat each pair @@ -98,7 +98,7 @@ MEAM::meam_force(int i, int eflag_global, int eflag_atom, int vflag_global, int eltj = fmap[type[j]]; scaleij = scale[type[i]][type[j]]; - printf("in j loop %d %d %d %d %f\n", i, j, eltj, type[j], scrfcn[fnoffset + jn]); + //printf("in j loop %d %d %d %d %f\n", i, j, eltj, type[j], scrfcn[fnoffset + jn]); if (!iszero(scrfcn[fnoffset + jn]) && eltj >= 0) { sij = scrfcn[fnoffset + jn] * fcpair[fnoffset + jn]; @@ -110,7 +110,7 @@ MEAM::meam_force(int i, int eflag_global, int eflag_atom, int vflag_global, int if (rij2 < this->cutforcesq) { rij = sqrt(rij2); recip = 1.0 / rij; - printf("forces, sij, rij = %f %f for i = %d and j = %d\n", sij, rij, i, j); + //printf("forces, sij, rij = %f %f for i = %d and j = %d\n", sij, rij, i, j); // Compute phi and phip ind = this->eltind[elti][eltj]; pp = rij * this->rdrar; @@ -120,12 +120,14 @@ MEAM::meam_force(int i, int eflag_global, int eflag_atom, int vflag_global, int pp = std::min(pp, 1.0); phi = ((this->phirar3[ind][kk] * pp + this->phirar2[ind][kk]) * pp + this->phirar1[ind][kk]) * pp + this->phirar[ind][kk]; phip = (this->phirar6[ind][kk] * pp + this->phirar5[ind][kk]) * pp + this->phirar4[ind][kk]; + /* printf("--- phi and phip i = %d ---\n", i); printf("phi %f %f\n", this->phirar3[ind][kk], this->phirar2[ind][kk]); printf("phi2 %f %f\n", this->phirar1[ind][kk], this->phirar[ind][kk]); printf("phip %f %f\n", this->phirar6[ind][kk], this->phirar5[ind][kk]); printf("phip2 %f\n", this->phirar4[ind][kk]); printf("-----"); + */ if (eflag_either != 0) { double phi_sc = phi * scaleij; @@ -212,9 +214,11 @@ MEAM::meam_force(int i, int eflag_global, int eflag_atom, int vflag_global, int const double t2mj = this->t2_meam[eltj]; const double t3mj = this->t3_meam[eltj]; - // ialloy mod not needed in msmeam according to fortran code - // if we remove this, need to multiply rhoali and rhoalj above by t vars - if (this->ialloy == 1) { + // ialloy mod not needed in msmeam, but similarity here is that we need to multiply rhos by t + // msmeam fortran code accomplishes this by multiplying rho's with t's above, like we did + // with rhoa1mj, rhoa2mj, etc + + if (this->ialloy == 1 || this->msmeamflag) { rhoa1j *= t1mj; rhoa2j *= t2mj; rhoa3j *= t3mj; @@ -229,10 +233,12 @@ MEAM::meam_force(int i, int eflag_global, int eflag_atom, int vflag_global, int drhoa3i *= t3mi; } + /* printf("forces, rhoali %f %f %f %f\n", rhoa0i, rhoa1i, rhoa2i, rhoa3i); printf("forces, rhoalmi %f %f %f\n", rhoa1mi, rhoa2mi, rhoa3mi); printf("forces, rhoalj %f %f %f %f\n", rhoa0j, rhoa1j, rhoa2j, rhoa3j); printf("forces, rhoalmj %f %f %f\n", rhoa1mj, rhoa2mj, rhoa3mj); + */ nv2 = 0; nv3 = 0; @@ -288,7 +294,7 @@ MEAM::meam_force(int i, int eflag_global, int eflag_atom, int vflag_global, int arg = delij[n] * delij[p] * this->v2D[nv2]; arg1i2m = arg1i2m + arho2m[i][nv2] * arg; arg1j2m = arg1j2m + arho2m[j][nv2] * arg; - printf(" arho2m[j][nv2] %d %d %f\n", j, nv2, arho2m[j][nv2]); + //printf(" arho2m[j][nv2] %d %d %f\n", j, nv2, arho2m[j][nv2]); nv2 = nv2 + 1; } arg1i1m = arg1i1m - arho1m[i][n] * delij[n]; @@ -298,10 +304,12 @@ MEAM::meam_force(int i, int eflag_global, int eflag_atom, int vflag_global, int } } + /* printf("rho1 terms arg1i %d %d %f %f\n", i, j, arg1i1, arg1i1m); printf("rho1 terms arg1j %d %d %f %f\n", i, j, arg1j1, arg1j1m); printf("rho1 terms rhoa1j %d %d %f %f\n", i, j, rhoa1j, rhoa1mj); printf("rho1 terms drhoa1j %d %d %f %f\n", i,j, drhoa1j, drhoa1mj); + */ // rho0 terms drho0dr1 = drhoa0j * sij; @@ -374,8 +382,8 @@ MEAM::meam_force(int i, int eflag_global, int eflag_atom, int vflag_global, int a2 = 2 * sij / rij2; drho2mdr1 = a2 * (drhoa2mj - 2 * rhoa2mj / rij) * arg1i2m - 2.0 / 3.0 * arho2mb[i] * drhoa2mj * sij; drho2mdr2 = a2 * (drhoa2mi - 2 * rhoa2mi / rij) * arg1j2m - 2.0 / 3.0 * arho2mb[j] * drhoa2mi * sij; - printf(" drho2mdr arho %d %d %f %f\n", i, j, arho2mb[i], arho2mb[j]); - printf(" drhoa rhoa arg1 %f %f %f\n", drhoa2mi, rhoa2mi, arg1j2m); + //printf(" drho2mdr arho %d %d %f %f\n", i, j, arho2mb[i], arho2mb[j]); + //printf(" drhoa rhoa arg1 %f %f %f\n", drhoa2mi, rhoa2mi, arg1j2m); a2 = 4 * sij / rij2; for (m = 0; m < 3; m++) { drho2mdrm1[m] = 0.0; @@ -388,13 +396,13 @@ MEAM::meam_force(int i, int eflag_global, int eflag_atom, int vflag_global, int drho2mdrm2[m] = -a2 * rhoa2mi * drho2mdrm2[m]; } - printf("forces drho2 %f %f %f %f\n", drho2drm1[0], drho2mdrm1[0], drho2drm2[0], drho2mdrm2[0]); + //printf("forces drho2 %f %f %f %f\n", drho2drm1[0], drho2mdrm1[0], drho2drm2[0], drho2mdrm2[0]); // rho3m terms rij3 = rij * rij2; a3 = 2 * sij / rij3; a3a = 6.0 / 5.0 * sij / rij; - printf("args %f %f %f %f\n", arg1i3m, arg3i3m, arg1j3m, arg3j3m); + //printf("args %f %f %f %f\n", arg1i3m, arg3i3m, arg1j3m, arg3j3m); drho3mdr1 = a3 * (drhoa3mj - 3 * rhoa3mj / rij) * arg1i3m - a3a * (drhoa3mj - rhoa3mj / rij) * arg3i3m; drho3mdr2 = a3 * (drhoa3mi - 3 * rhoa3mi / rij) * arg1j3m - a3a * (drhoa3mi - rhoa3mi / rij) * arg3j3m; drho3mdr1 *= -1.0; @@ -402,7 +410,7 @@ MEAM::meam_force(int i, int eflag_global, int eflag_atom, int vflag_global, int a3 = 6 * sij / rij3; a3a = 6 * sij / (5 * rij); - printf("forces drho3mdr %f %f\n", drho3mdr1, drho3mdr2); + //printf("forces drho3mdr %f %f\n", drho3mdr1, drho3mdr2); for (m = 0; m < 3; m++) { drho3mdrm1[m] = 0.0; drho3mdrm2[m] = 0.0; @@ -426,11 +434,12 @@ MEAM::meam_force(int i, int eflag_global, int eflag_atom, int vflag_global, int drho3mdrm1[m] = 0.0; drho3mdrm2[m] = 0.0; } - + /* printf("rho2 terms drho2dr1 %f %f\n", drho2dr1, drho2mdr1); printf("rho2 terms drho2dr2 %f %f\n", drho2dr2, drho2mdr2); printf("rho2 terms %f %f\n", drho2drm1[0], drho2mdrm1[0]); printf("rho2 terms %f %f\n", drho2drm2[0], drho2mdrm2[0]); + */ // compute derivatives of weighting functions t wrt rij @@ -509,6 +518,7 @@ MEAM::meam_force(int i, int eflag_global, int eflag_atom, int vflag_global, int get_shpfcn(this->lattce_meam[eltj][eltj], this->stheta_meam[elti][elti], this->ctheta_meam[elti][elti], shpj); if (this->msmeamflag){ + /* printf("test rhodg %f %f %f\n", dgamma1[i], dgamma2[i], dgamma3[i]); printf("test rhod dtdr %f %f %f\n", dt1dr1, dt2dr1, dt3dr1); printf("test rhod rho %f %f %f\n", rho1[i], rho2[i], rho3[i]); @@ -520,6 +530,7 @@ MEAM::meam_force(int i, int eflag_global, int eflag_atom, int vflag_global, int printf(" %f %f %f\n", dgamma1[i], dgamma2[i], dgamma3[i]); printf(" %f %f %f\n", drho1dr2, drho2dr2, drho3dr2); printf(" %f %f %f\n", drho1mdr2, drho2mdr2, drho3mdr2); + */ drhodr1 = dgamma1[i] * drho0dr1 + dgamma2[i] * (dt1dr1 * rho1[i] + t1i * (drho1dr1 - drho1mdr1) + dt2dr1 * rho2[i] + t2i * (drho2dr1 - drho2mdr1) + @@ -530,7 +541,7 @@ MEAM::meam_force(int i, int eflag_global, int eflag_atom, int vflag_global, int dt2dr2 * rho2[j] + t2j * (drho2dr2 - drho2mdr2) + dt3dr2 * rho3[j] + t3j * (drho3dr2 - drho3mdr2)) - dgamma3[j] * (shpj[0] * dt1dr2 + shpj[1] * dt2dr2 + shpj[2] * dt3dr2); - printf(" drhodr1&2 %f %f\n", drhodr1, drhodr2); + //printf(" drhodr1&2 %f %f\n", drhodr1, drhodr2); for (m = 0; m < 3; m++) { drhodrm1[m] = 0.0; drhodrm2[m] = 0.0; @@ -559,12 +570,12 @@ MEAM::meam_force(int i, int eflag_global, int eflag_atom, int vflag_global, int } } - printf("rhod terms drhodr1 %f %f\n", drhodr1, drhodr2); - printf("rhod terms drhodr1 %f %f\n", drhodrm1[0], drhodrm2[0]); + //printf("rhod terms drhodr1 %f %f\n", drhodr1, drhodr2); + //printf("rhod terms drhodr1 %f %f\n", drhodrm1[0], drhodrm2[0]); // Compute derivatives wrt sij, but only if necessary if (!iszero(dscrfcn[fnoffset + jn])) { - printf("computing derivatives wrt sij\n"); + //printf("computing derivatives wrt sij\n"); drho0ds1 = rhoa0j; drho0ds2 = rhoa0i; a1 = 2.0 / rij; @@ -598,6 +609,7 @@ MEAM::meam_force(int i, int eflag_global, int eflag_atom, int vflag_global, int } if (this->ialloy == 1) { + a1i = fdiv_zero(rhoa0j, tsq_ave[i][0]); a1j = fdiv_zero(rhoa0i, tsq_ave[j][0]); a2i = fdiv_zero(rhoa0j, tsq_ave[i][1]); @@ -667,11 +679,13 @@ MEAM::meam_force(int i, int eflag_global, int eflag_atom, int vflag_global, int dUdrij = phip * sij + frhop[i] * drhodr1 + frhop[j] * drhodr2; dUdsij = 0.0; + printf("--- compute derivatives of energy wrt ...\n"); printf("phip sij %f %f\n", phip,sij); printf("drhodr %f %f\n", drhodr1, drhodr2); printf("fp %f %f\n", frhop[i], frhop[j]); printf("dUdrij %f\n", dUdrij); + if (!iszero(dscrfcn[fnoffset + jn])) { dUdsij = phi + frhop[i] * drhods1 + frhop[j] * drhods2; } diff --git a/src/MEAM/meam_setup_done.cpp b/src/MEAM/meam_setup_done.cpp index c3dfc8c7fc..0335b7ae3d 100644 --- a/src/MEAM/meam_setup_done.cpp +++ b/src/MEAM/meam_setup_done.cpp @@ -377,7 +377,7 @@ double MEAM::phi_meam(double r, int a, int b) // calculate average weighting factors for the reference structure if (this->lattce_meam[a][b] == C11) { - // note ialloy and t_ave not used in MS-MEAM + // note ialloy and t_ave not used in msmeam if (this->ialloy == 2) { t11av = this->t1_meam[a]; t12av = this->t1_meam[b]; @@ -399,7 +399,7 @@ double MEAM::phi_meam(double r, int a, int b) get_tavref(&t11av, &t21av, &t31av, &t12av, &t22av, &t32av, this->t1_meam[a], this->t2_meam[a], this->t3_meam[a], this->t1_meam[b], this->t2_meam[b], this->t3_meam[b], r, a, b, this->lattce_meam[a][b]); - // msmeam - call twice with different sets of variables + // with msmeam call twice with different sets of variables if (this->msmeamflag){ get_tavref(&t1m1av, &t2m1av, &t3m1av, &t1m2av, &t2m2av, &t3m2av, this->t1m_meam[a], this->t2m_meam[a], this->t3m_meam[a], this->t1m_meam[b], this->t2m_meam[b], this->t3m_meam[b], r, a, b, @@ -731,6 +731,7 @@ void MEAM::get_densref(double r, int a, int b, double* rho01, double* rho11, dou rhoa22 = this->rho0_meam[b] * MathSpecial::fm_exp(-this->beta2_meam[b] * a2); rhoa32 = this->rho0_meam[b] * MathSpecial::fm_exp(-this->beta3_meam[b] * a2); */ + /* rhoa11 = this->rho0_meam[a] * this->t1_meam[a] * MathSpecial::fm_exp(-this->beta1_meam[a] * a1); rhoa21 = this->rho0_meam[a] * this->t2_meam[a] * MathSpecial::fm_exp(-this->beta2_meam[a] * a1); rhoa31 = this->rho0_meam[a] * this->t3_meam[a] * MathSpecial::fm_exp(-this->beta3_meam[a] * a1); @@ -738,7 +739,17 @@ void MEAM::get_densref(double r, int a, int b, double* rho01, double* rho11, dou rhoa12 = this->rho0_meam[b] * this->t1_meam[b] * MathSpecial::fm_exp(-this->beta1_meam[b] * a2); rhoa22 = this->rho0_meam[b] * this->t2_meam[b] * MathSpecial::fm_exp(-this->beta2_meam[b] * a2); rhoa32 = this->rho0_meam[b] * this->t3_meam[b] * MathSpecial::fm_exp(-this->beta3_meam[b] * a2); + */ if (this->msmeamflag){ + // in msmeam, the rho variables are multiplied by t here since we don't use ialloy + rhoa11 = this->rho0_meam[a] * this->t1_meam[a] * MathSpecial::fm_exp(-this->beta1_meam[a] * a1); + rhoa21 = this->rho0_meam[a] * this->t2_meam[a] * MathSpecial::fm_exp(-this->beta2_meam[a] * a1); + rhoa31 = this->rho0_meam[a] * this->t3_meam[a] * MathSpecial::fm_exp(-this->beta3_meam[a] * a1); + rhoa02 = this->rho0_meam[b] * MathSpecial::fm_exp(-this->beta0_meam[b] * a2); + rhoa12 = this->rho0_meam[b] * this->t1_meam[b] * MathSpecial::fm_exp(-this->beta1_meam[b] * a2); + rhoa22 = this->rho0_meam[b] * this->t2_meam[b] * MathSpecial::fm_exp(-this->beta2_meam[b] * a2); + rhoa32 = this->rho0_meam[b] * this->t3_meam[b] * MathSpecial::fm_exp(-this->beta3_meam[b] * a2); + // msmeam specific rho vars rhoa1m1 = this->rho0_meam[a] * this->t1m_meam[a] * MathSpecial::fm_exp(-this->beta1m_meam[a] * a1); rhoa2m1 = this->rho0_meam[a] * this->t2m_meam[a] * MathSpecial::fm_exp(-this->beta2m_meam[a] * a1); rhoa3m1 = this->rho0_meam[a] * this->t3m_meam[a] * MathSpecial::fm_exp(-this->beta3m_meam[a] * a1); @@ -746,6 +757,15 @@ void MEAM::get_densref(double r, int a, int b, double* rho01, double* rho11, dou rhoa2m2 = this->rho0_meam[b] * this->t2m_meam[b] * MathSpecial::fm_exp(-this->beta2m_meam[b] * a2); rhoa3m2 = this->rho0_meam[b] * this->t3m_meam[b] * MathSpecial::fm_exp(-this->beta3m_meam[b] * a2); } + else{ + rhoa11 = this->rho0_meam[a] * MathSpecial::fm_exp(-this->beta1_meam[a] * a1); + rhoa21 = this->rho0_meam[a] * MathSpecial::fm_exp(-this->beta2_meam[a] * a1); + rhoa31 = this->rho0_meam[a] * MathSpecial::fm_exp(-this->beta3_meam[a] * a1); + rhoa02 = this->rho0_meam[b] * MathSpecial::fm_exp(-this->beta0_meam[b] * a2); + rhoa12 = this->rho0_meam[b] * MathSpecial::fm_exp(-this->beta1_meam[b] * a2); + rhoa22 = this->rho0_meam[b] * MathSpecial::fm_exp(-this->beta2_meam[b] * a2); + rhoa32 = this->rho0_meam[b] * MathSpecial::fm_exp(-this->beta3_meam[b] * a2); + } //printf("rhoa0 %f %f\n", rhoa01, rhoa02); //printf("rhoai1 %f %f %f\n", rhoa11, rhoa21, rhoa31); //printf("rhoai2 %f %f %f\n", rhoa12, rhoa22, rhoa32); @@ -842,9 +862,11 @@ void MEAM::get_densref(double r, int a, int b, double* rho01, double* rho11, dou case L12: *rho01 = 8 * rhoa01 + 4 * rhoa02; *rho02 = 12 * rhoa01; - // not modified by choice with msmeam; use of t's inconsistent with above + // don't use with msmeam; use of t's inconsistent with above + // also we should not use ialloy with msmeam + // possibly make ialloy not an option when using msmeam if (this->ialloy == 1 && !this->msmeamflag) { - printf("^^^ ialloy=1 in getdensref\n"); + //printf("^^^ ialloy=1 in getdensref\n"); *rho21 = 8. / 3. * MathSpecial::square(rhoa21 * this->t2_meam[a] - rhoa22 * this->t2_meam[b]); denom = 8 * rhoa01 * MathSpecial::square(this->t2_meam[a]) + 4 * rhoa02 * MathSpecial::square(this->t2_meam[b]); if (denom > 0.) diff --git a/src/MEAM/meam_setup_global.cpp b/src/MEAM/meam_setup_global.cpp index 0ff64b06bc..859de87d5a 100644 --- a/src/MEAM/meam_setup_global.cpp +++ b/src/MEAM/meam_setup_global.cpp @@ -40,11 +40,6 @@ MEAM::meam_setup_global(int nelt, lattice_t* lat, int* ielement, double* /*atwt* double *t3m) { - printf("^^^^^ meam setup global msmeamflag: %d\n", this->msmeamflag); - if (this->msmeamflag){ - printf("b1m[0]: %f\n", b1m[0]); - } - int i; double tmplat[maxelt]; @@ -60,7 +55,7 @@ MEAM::meam_setup_global(int nelt, lattice_t* lat, int* ielement, double* /*atwt* this->beta2_meam[i] = b2[i]; this->beta3_meam[i] = b3[i]; if (this->msmeamflag){ - printf("^^^^^ Filling betam\n"); + //printf("^^^^^ Filling betam\n"); this->beta1m_meam[i] = b1m[i]; this->beta2m_meam[i] = b2m[i]; this->beta3m_meam[i] = b3m[i]; diff --git a/src/MEAM/pair_meam.cpp b/src/MEAM/pair_meam.cpp index a53b0e3cf8..0aec140eed 100644 --- a/src/MEAM/pair_meam.cpp +++ b/src/MEAM/pair_meam.cpp @@ -134,7 +134,7 @@ void PairMEAM::compute(int eflag, int vflag) int offset = 0; errorflag = 0; - printf("before meam_dens_init %f %f\n", meam_inst->arho2m[0][0], meam_inst->arho3m[0][0]); + //printf("before meam_dens_init %f %f\n", meam_inst->arho2m[0][0], meam_inst->arho3m[0][0]); for (ii = 0; ii < inum_half; ii++) { i = ilist_half[ii]; @@ -144,18 +144,18 @@ void PairMEAM::compute(int eflag, int vflag) offset); offset += numneigh_half[i]; } - printf("after meam_dens_init %f %f\n", meam_inst->arho2m[0][0], meam_inst->arho3m[0][0]); - printf("before revcomm %f %f\n", meam_inst->arho2m[0][0], meam_inst->arho3m[0][0]); + //printf("after meam_dens_init %f %f\n", meam_inst->arho2m[0][0], meam_inst->arho3m[0][0]); + //printf("before revcomm %f %f\n", meam_inst->arho2m[0][0], meam_inst->arho3m[0][0]); comm->reverse_comm(this); - printf("after revcomm %f %f\n", meam_inst->arho2m[0][0], meam_inst->arho3m[0][0]); + //printf("after revcomm %f %f\n", meam_inst->arho2m[0][0], meam_inst->arho3m[0][0]); meam_inst->meam_dens_final(nlocal,eflag_either,eflag_global,eflag_atom, &eng_vdwl,eatom,ntype,type,map,scale,errorflag); if (errorflag) error->one(FLERR,"MEAM library error {}",errorflag); - printf("before forcomm %f %f\n", meam_inst->arho2m[0][0], meam_inst->arho3m[0][0]); + //printf("before forcomm %f %f\n", meam_inst->arho2m[0][0], meam_inst->arho3m[0][0]); comm->forward_comm(this); - printf("after forcomm %f %f\n", meam_inst->arho2m[0][0], meam_inst->arho3m[0][0]); + //printf("after forcomm %f %f\n", meam_inst->arho2m[0][0], meam_inst->arho3m[0][0]); offset = 0;