diff --git a/doc/src/pair_born.txt b/doc/src/pair_born.txt index 398b9b1717..caf901485d 100644 --- a/doc/src/pair_born.txt +++ b/doc/src/pair_born.txt @@ -19,6 +19,8 @@ pair_style born/coul/msm/omp command :h3 pair_style born/coul/wolf command :h3 pair_style born/coul/wolf/gpu command :h3 pair_style born/coul/wolf/omp command :h3 +pair_style born/coul/dsf command :h3 +pair_style born/coul/dsf/cs command :h3 [Syntax:] @@ -37,7 +39,11 @@ args = list of arguments for a particular style :ul {born/coul/wolf} args = alpha cutoff (cutoff2) alpha = damping parameter (inverse distance units) cutoff = global cutoff for non-Coulombic (and Coulombic if only 1 arg) (distance units) - cutoff2 = global cutoff for Coulombic (optional) (distance units) :pre + cutoff2 = global cutoff for Coulombic (optional) (distance units) + {born/coul/dsf} or {born/coul/dsf/cs} args = alpha cutoff (cutoff2) + alpha = damping parameter (inverse distance units) + cutoff = global cutoff for non-Coulombic (and Coulombic if only 1 arg) (distance units) + cutoff2 = global cutoff for Coulombic (distance units) :pre [Examples:] @@ -62,6 +68,11 @@ pair_style born/coul/wolf 0.25 10.0 9.0 pair_coeff * * 6.08 0.317 2.340 24.18 11.51 pair_coeff 1 1 6.08 0.317 2.340 24.18 11.51 :pre +pair_style born/coul/dsf 0.1 10.0 12.0 +pair_coeff * * 0.0 1.00 0.00 0.00 0.00 +pair_coeff 1 1 480.0 0.25 0.00 1.05 0.50 :pre + + [Description:] The {born} style computes the Born-Mayer-Huggins or Tosi/Fumi @@ -90,10 +101,15 @@ term. The {born/coul/wolf} style adds a Coulombic term as described for the Wolf potential in the "coul/wolf"_pair_coul.html pair style. +The {born/coul/dsf} style computes the Coulomb contribution +with the damped shifted force model as in the +"coul/dsf"_pair_coul.html style. + Style {born/coul/long/cs} is identical to {born/coul/long} except that a term is added for the "core/shell model"_Section_howto.html#howto_25 to allow charges on core and shell particles to be separated by r = -0.0. +0.0. The same correction is introduced for {born/coul/dsf/cs} style +which is identical to {born/coul/dsf}. Note that these potentials are related to the "Buckingham potential"_pair_buck.html. @@ -116,9 +132,10 @@ The second coefficient, rho, must be greater than zero. The last coefficient is optional. If not specified, the global A,C,D cutoff specified in the pair_style command is used. -For {born/coul/long} and {born/coul/wolf} no Coulombic cutoff can be -specified for an individual I,J type pair. All type pairs use the -same global Coulombic cutoff specified in the pair_style command. +For {born/coul/long}, {born/coul/wolf} and {born/coul/dsf} no +Coulombic cutoff can be specified for an individual I,J type pair. +All type pairs use the same global Coulombic cutoff specified in +the pair_style command. :line diff --git a/doc/src/pair_cs.txt b/doc/src/pair_cs.txt index b5fcfdc39c..3e1390049a 100644 --- a/doc/src/pair_cs.txt +++ b/doc/src/pair_cs.txt @@ -8,19 +8,24 @@ pair_style born/coul/long/cs command :h3 pair_style buck/coul/long/cs command :h3 +pair_style born/coul/dsf/cs command :h3 [Syntax:] pair_style style args :pre -style = {born/coul/long/cs} or {buck/coul/long/cs} +style = {born/coul/long/cs} or {buck/coul/long/cs} or {born/coul/dsf/cs} args = list of arguments for a particular style :ul {born/coul/long/cs} args = cutoff (cutoff2) cutoff = global cutoff for non-Coulombic (and Coulombic if only 1 arg) (distance units) cutoff2 = global cutoff for Coulombic (optional) (distance units) {buck/coul/long/cs} args = cutoff (cutoff2) cutoff = global cutoff for Buckingham (and Coulombic if only 1 arg) (distance units) - cutoff2 = global cutoff for Coulombic (optional) (distance units) :pre + cutoff2 = global cutoff for Coulombic (optional) (distance units) + {born/coul/dsf/cs} args = alpha cutoff (cutoff2) + alpha = damping parameter (inverse distance units) + cutoff = global cutoff for non-Coulombic (and Coulombic if only 1 arg) (distance units) + cutoff2 = global cutoff for Coulombic (distance units) :pre [Examples:] @@ -32,6 +37,10 @@ pair_style buck/coul/long/cs 10.0 8.0 pair_coeff * * 100.0 1.5 200.0 pair_coeff 1 1 100.0 1.5 200.0 9.0 :pre +pair_style born/coul/dsf/cs 0.1 10.0 12.0 +pair_coeff * * 0.0 1.00 0.00 0.00 0.00 +pair_coeff 1 1 480.0 0.25 0.00 1.05 0.50 :pre + [Description:] These pair styles are designed to be used with the adiabatic @@ -39,7 +48,7 @@ core/shell model of "(Mitchell and Finchham)"_#MitchellFinchham. See "Section 6.25"_Section_howto.html#howto_25 of the manual for an overview of the model as implemented in LAMMPS. -These pair styles are identical to the "pair_style +The styles with a {coul/long} term are identical to the "pair_style born/coul/long"_pair_born.html and "pair_style buck/coul/long"_pair_buck.html styles, except they correctly treat the special case where the distance between two charged core and shell @@ -63,6 +72,14 @@ where C is an energy-conversion constant, Qi and Qj are the charges on the core and shell, epsilon is the dielectric constant and r_min is the minimal distance. +The pair style {born/coul/dsf/cs} is identical to the +"pair_style born/coul/dsf"_pair_born.html style, which uses the +the damped shifted force model as in "coul/dsf"_pair_coul.html +to compute the Coulomb contribution. This approach does not require +a long-range solver, thus the only correction is the addition of a +minimal distance to avoid the possible r = 0.0 case for a +core/shell pair. + [Restrictions:] These pair styles are part of the CORESHELL package. They are only diff --git a/examples/coreshell/in.coreshell.dsf b/examples/coreshell/in.coreshell.dsf new file mode 100644 index 0000000000..5c9d6c7c03 --- /dev/null +++ b/examples/coreshell/in.coreshell.dsf @@ -0,0 +1,71 @@ +# Testsystem for core-shell model compared to Mitchel and Finchham +# Hendrik Heenen, June 2014 + +# ------------------------ INITIALIZATION ---------------------------- + +units metal +dimension 3 +boundary p p p +atom_style full + +# ----------------------- ATOM DEFINITION ---------------------------- + +fix csinfo all property/atom i_CSID +read_data data.coreshell fix csinfo NULL CS-Info + +group cores type 1 2 +group shells type 3 4 + +neighbor 2.0 bin +comm_modify vel yes + +# ------------------------ FORCE FIELDS ------------------------------ + +pair_style born/coul/dsf/cs 0.1 20.0 20.0 # A, rho, sigma=0, C, D +pair_coeff * * 0.0 1.000 0.00 0.00 0.00 +pair_coeff 3 3 487.0 0.23768 0.00 1.05 0.50 #Na-Na +pair_coeff 3 4 145134.0 0.23768 0.00 6.99 8.70 #Na-Cl +pair_coeff 4 4 405774.0 0.23768 0.00 72.40 145.40 #Cl-Cl + +bond_style harmonic +bond_coeff 1 63.014 0.0 +bond_coeff 2 25.724 0.0 + +# ------------------------ Equilibration Run ------------------------------- + +reset_timestep 0 + +thermo 50 +thermo_style custom step etotal pe ke temp press & + epair evdwl ecoul elong ebond fnorm fmax vol + +compute CSequ all temp/cs cores shells + +# output via chunk method + +#compute prop all property/atom i_CSID +#compute cs_chunk all chunk/atom c_prop +#compute cstherm all temp/chunk cs_chunk temp internal com yes cdof 3.0 +#fix ave_chunk all ave/time 100 1 100 c_cstherm file chunk.dump mode vector + +thermo_modify temp CSequ + +# velocity bias option + +velocity all create 1427 134 dist gaussian mom yes rot no bias yes temp CSequ +velocity all scale 1427 temp CSequ + +fix thermoberendsen all temp/berendsen 1427 1427 0.4 +fix nve all nve +fix_modify thermoberendsen temp CSequ + +# 2 fmsec timestep + +timestep 0.002 +run 500 + +unfix thermoberendsen + +# ------------------------ Dynamic Run ------------------------------- + +run 1000 diff --git a/examples/coreshell/log.8Nov16.coreshell.dsf.1 b/examples/coreshell/log.8Nov16.coreshell.dsf.1 new file mode 100644 index 0000000000..f58970fff0 --- /dev/null +++ b/examples/coreshell/log.8Nov16.coreshell.dsf.1 @@ -0,0 +1,185 @@ +LAMMPS (27 Oct 2016) +# Testsystem for core-shell model compared to Mitchel and Finchham +# Hendrik Heenen, June 2014 + +# ------------------------ INITIALIZATION ---------------------------- + +units metal +dimension 3 +boundary p p p +atom_style full + +# ----------------------- ATOM DEFINITION ---------------------------- + +fix csinfo all property/atom i_CSID +read_data data.coreshell fix csinfo NULL CS-Info + orthogonal box = (0 0 0) to (24.096 24.096 24.096) + 1 by 1 by 1 MPI processor grid + reading atoms ... + 432 atoms + scanning bonds ... + 1 = max bonds/atom + reading bonds ... + 216 bonds + 1 = max # of 1-2 neighbors + 0 = max # of 1-3 neighbors + 0 = max # of 1-4 neighbors + 1 = max # of special neighbors + +group cores type 1 2 +216 atoms in group cores +group shells type 3 4 +216 atoms in group shells + +neighbor 2.0 bin +comm_modify vel yes + +# ------------------------ FORCE FIELDS ------------------------------ + +pair_style born/coul/dsf/cs 0.1 20.0 20.0 # A, rho, sigma=0, C, D +pair_coeff * * 0.0 1.000 0.00 0.00 0.00 +pair_coeff 3 3 487.0 0.23768 0.00 1.05 0.50 #Na-Na +pair_coeff 3 4 145134.0 0.23768 0.00 6.99 8.70 #Na-Cl +pair_coeff 4 4 405774.0 0.23768 0.00 72.40 145.40 #Cl-Cl + +bond_style harmonic +bond_coeff 1 63.014 0.0 +bond_coeff 2 25.724 0.0 + +# ------------------------ Equilibration Run ------------------------------- + +reset_timestep 0 + +thermo 50 +thermo_style custom step etotal pe ke temp press epair evdwl ecoul elong ebond fnorm fmax vol + +compute CSequ all temp/cs cores shells + +# output via chunk method + +#compute prop all property/atom i_CSID +#compute cs_chunk all chunk/atom c_prop +#compute cstherm all temp/chunk cs_chunk temp internal com yes cdof 3.0 +#fix ave_chunk all ave/time 100 1 100 c_cstherm file chunk.dump mode vector + +thermo_modify temp CSequ + +# velocity bias option + +velocity all create 1427 134 dist gaussian mom yes rot no bias yes temp CSequ +Neighbor list info ... + 1 neighbor list requests + update every 1 steps, delay 10 steps, check yes + max neighbors/atom: 2000, page size: 100000 + master list distance cutoff = 22 + ghost atom cutoff = 22 + binsize = 11 -> bins = 3 3 3 +velocity all scale 1427 temp CSequ + +fix thermoberendsen all temp/berendsen 1427 1427 0.4 +fix nve all nve +fix_modify thermoberendsen temp CSequ + +# 2 fmsec timestep + +timestep 0.002 +run 500 +Memory usage per processor = 7.04355 Mbytes +Step TotEng PotEng KinEng Temp Press E_pair E_vdwl E_coul E_long E_bond Fnorm Fmax Volume + 0 -635.80596 -675.46362 39.657659 1427 -21302.622 -675.46362 1.6320365 -677.09565 0 0 1.3517686e-14 2.942091e-15 13990.5 + 50 -633.9898 -666.02679 32.03699 1152.7858 -4578.5681 -668.50431 37.800204 -706.30452 0 2.4775226 14.568073 4.3012389 13990.5 + 100 -631.89604 -661.96148 30.065442 1081.8436 -3536.6738 -664.61798 39.18583 -703.80381 0 2.6564973 14.677968 3.9051029 13990.5 + 150 -630.08723 -662.95879 32.871559 1182.816 -109.19506 -665.76772 46.247821 -712.01554 0 2.8089226 15.270039 2.9328953 13990.5 + 200 -628.55895 -663.97376 35.414806 1274.3296 -1748.35 -666.58439 41.738552 -708.32294 0 2.6106349 14.148282 3.1047826 13990.5 + 250 -627.28761 -661.92274 34.635123 1246.2743 -1280.4899 -664.917 43.045475 -707.96247 0 2.9942594 14.248617 2.4694705 13990.5 + 300 -626.6163 -663.65651 37.040209 1332.8164 -1887.9043 -666.35215 40.84964 -707.20179 0 2.6956373 13.142643 1.9263242 13990.5 + 350 -625.76781 -664.66441 38.896607 1399.6151 -1839.482 -667.47659 40.999206 -708.47579 0 2.8121749 13.601238 1.9262698 13990.5 + 400 -625.02586 -661.46042 36.434568 1311.0236 -868.2031 -664.40231 43.21398 -707.61629 0 2.9418875 14.945389 2.7493413 13990.5 + 450 -624.3278 -660.50844 36.180639 1301.8865 -2203.3944 -663.49896 40.008669 -703.50763 0 2.9905179 14.158866 1.7299899 13990.5 + 500 -623.56254 -661.33839 37.775849 1359.2869 -810.50736 -664.11652 42.993999 -707.11052 0 2.7781274 13.68709 2.9115277 13990.5 +Loop time of 10.7162 on 1 procs for 500 steps with 432 atoms + +Performance: 8.063 ns/day, 2.977 hours/ns, 46.658 timesteps/s +99.9% CPU use with 1 MPI tasks x no OpenMP threads + +MPI task timing breakdown: +Section | min time | avg time | max time |%varavg| %total +--------------------------------------------------------------- +Pair | 10.478 | 10.478 | 10.478 | 0.0 | 97.78 +Bond | 0.0029511 | 0.0029511 | 0.0029511 | 0.0 | 0.03 +Neigh | 0.14159 | 0.14159 | 0.14159 | 0.0 | 1.32 +Comm | 0.074382 | 0.074382 | 0.074382 | 0.0 | 0.69 +Output | 0.00054097 | 0.00054097 | 0.00054097 | 0.0 | 0.01 +Modify | 0.010588 | 0.010588 | 0.010588 | 0.0 | 0.10 +Other | | 0.007748 | | | 0.07 + +Nlocal: 432 ave 432 max 432 min +Histogram: 1 0 0 0 0 0 0 0 0 0 +Nghost: 9280 ave 9280 max 9280 min +Histogram: 1 0 0 0 0 0 0 0 0 0 +Neighs: 297636 ave 297636 max 297636 min +Histogram: 1 0 0 0 0 0 0 0 0 0 + +Total # of neighbors = 297636 +Ave neighs/atom = 688.972 +Ave special neighs/atom = 1 +Neighbor list builds = 20 +Dangerous builds = 0 + +unfix thermoberendsen + +# ------------------------ Dynamic Run ------------------------------- + +run 1000 +Memory usage per processor = 7.04355 Mbytes +Step TotEng PotEng KinEng Temp Press E_pair E_vdwl E_coul E_long E_bond Fnorm Fmax Volume + 500 -623.56254 -661.33839 37.775849 1359.2869 -810.50736 -664.11652 42.993999 -707.11052 0 2.7781274 13.68709 2.9115277 13990.5 + 550 -623.5004 -660.74472 37.244326 1340.1611 -1413.4326 -663.99669 41.875014 -705.8717 0 3.2519651 15.097948 2.278405 13990.5 + 600 -623.46963 -659.61729 36.147655 1300.6997 -521.50578 -662.54994 43.956071 -706.50601 0 2.9326492 14.99649 2.6334959 13990.5 + 650 -623.49291 -661.50698 38.014069 1367.8588 -1230.0925 -664.21074 42.027844 -706.23859 0 2.7037578 13.982308 1.6247207 13990.5 + 700 -623.4913 -660.11564 36.62434 1317.8522 -727.89052 -663.24921 43.413397 -706.66261 0 3.1335699 15.009937 2.0563966 13990.5 + 750 -623.50292 -657.95982 34.4569 1239.8613 636.46644 -661.16971 46.539267 -707.70898 0 3.2098934 15.25993 2.1864622 13990.5 + 800 -623.5176 -659.92032 36.402711 1309.8773 -912.75799 -662.84989 42.668309 -705.5182 0 2.9295708 13.577516 2.0006099 13990.5 + 850 -623.44098 -660.92727 37.486295 1348.8679 -550.40358 -664.08308 43.667245 -707.75033 0 3.1558098 14.836208 2.279198 13990.5 + 900 -623.46361 -661.21737 37.753765 1358.4923 1267.8647 -664.52195 47.67284 -712.19479 0 3.3045765 15.058502 1.886141 13990.5 + 950 -623.50114 -660.58464 37.083492 1334.3739 1754.7359 -663.48186 48.70363 -712.18549 0 2.897226 15.519042 2.2654928 13990.5 + 1000 -623.50161 -660.02915 36.527539 1314.369 228.76104 -663.31152 45.374099 -708.68562 0 3.2823685 14.783709 2.4201134 13990.5 + 1050 -623.45985 -660.57417 37.114321 1335.4832 -1490.604 -663.75391 41.258878 -705.01279 0 3.1797391 14.250262 2.3153255 13990.5 + 1100 -623.51051 -661.20338 37.692871 1356.3011 1791.7899 -664.01042 48.626451 -712.63687 0 2.807039 15.559872 3.184101 13990.5 + 1150 -623.51067 -663.19545 39.684776 1427.9758 1023.0584 -666.07723 46.5628 -712.64003 0 2.8817804 13.895322 2.3950292 13990.5 + 1200 -623.49625 -659.6715 36.175253 1301.6927 1600.2805 -662.62259 48.522365 -711.14495 0 2.9510854 15.567834 2.1677651 13990.5 + 1250 -623.48282 -660.56735 37.084533 1334.4113 -871.67341 -663.86673 42.560699 -706.42743 0 3.2993759 14.569539 2.0093709 13990.5 + 1300 -623.47744 -663.63125 40.153811 1444.853 1343.7147 -666.39564 47.104842 -713.50048 0 2.7643857 14.186019 1.4599359 13990.5 + 1350 -623.49121 -661.42731 37.936096 1365.0531 589.73669 -664.46099 45.947687 -710.40867 0 3.0336821 14.801223 2.7486556 13990.5 + 1400 -623.50803 -660.03912 36.53109 1314.4968 362.97431 -663.24909 45.772904 -709.02199 0 3.2099708 14.566488 1.9170714 13990.5 + 1450 -623.51243 -659.65548 36.143052 1300.534 2853.0755 -663.0534 51.355353 -714.40875 0 3.3979157 15.890282 2.5251359 13990.5 + 1500 -623.51621 -661.87741 38.361201 1380.3496 740.04973 -665.00896 46.208742 -711.2177 0 3.1315492 15.168927 2.4710846 13990.5 +Loop time of 22.2766 on 1 procs for 1000 steps with 432 atoms + +Performance: 7.757 ns/day, 3.094 hours/ns, 44.890 timesteps/s +99.9% CPU use with 1 MPI tasks x no OpenMP threads + +MPI task timing breakdown: +Section | min time | avg time | max time |%varavg| %total +--------------------------------------------------------------- +Pair | 21.8 | 21.8 | 21.8 | 0.0 | 97.86 +Bond | 0.005852 | 0.005852 | 0.005852 | 0.0 | 0.03 +Neigh | 0.30423 | 0.30423 | 0.30423 | 0.0 | 1.37 +Comm | 0.14388 | 0.14388 | 0.14388 | 0.0 | 0.65 +Output | 0.0010855 | 0.0010855 | 0.0010855 | 0.0 | 0.00 +Modify | 0.0064189 | 0.0064189 | 0.0064189 | 0.0 | 0.03 +Other | | 0.01527 | | | 0.07 + +Nlocal: 432 ave 432 max 432 min +Histogram: 1 0 0 0 0 0 0 0 0 0 +Nghost: 9318 ave 9318 max 9318 min +Histogram: 1 0 0 0 0 0 0 0 0 0 +Neighs: 297131 ave 297131 max 297131 min +Histogram: 1 0 0 0 0 0 0 0 0 0 + +Total # of neighbors = 297131 +Ave neighs/atom = 687.803 +Ave special neighs/atom = 1 +Neighbor list builds = 44 +Dangerous builds = 0 +Total wall time: 0:00:33 diff --git a/examples/coreshell/log.8Nov16.coreshell.dsf.2 b/examples/coreshell/log.8Nov16.coreshell.dsf.2 new file mode 100644 index 0000000000..94da2f5832 --- /dev/null +++ b/examples/coreshell/log.8Nov16.coreshell.dsf.2 @@ -0,0 +1,185 @@ +LAMMPS (27 Oct 2016) +# Testsystem for core-shell model compared to Mitchel and Finchham +# Hendrik Heenen, June 2014 + +# ------------------------ INITIALIZATION ---------------------------- + +units metal +dimension 3 +boundary p p p +atom_style full + +# ----------------------- ATOM DEFINITION ---------------------------- + +fix csinfo all property/atom i_CSID +read_data data.coreshell fix csinfo NULL CS-Info + orthogonal box = (0 0 0) to (24.096 24.096 24.096) + 1 by 1 by 2 MPI processor grid + reading atoms ... + 432 atoms + scanning bonds ... + 1 = max bonds/atom + reading bonds ... + 216 bonds + 1 = max # of 1-2 neighbors + 0 = max # of 1-3 neighbors + 0 = max # of 1-4 neighbors + 1 = max # of special neighbors + +group cores type 1 2 +216 atoms in group cores +group shells type 3 4 +216 atoms in group shells + +neighbor 2.0 bin +comm_modify vel yes + +# ------------------------ FORCE FIELDS ------------------------------ + +pair_style born/coul/dsf/cs 0.1 20.0 20.0 # A, rho, sigma=0, C, D +pair_coeff * * 0.0 1.000 0.00 0.00 0.00 +pair_coeff 3 3 487.0 0.23768 0.00 1.05 0.50 #Na-Na +pair_coeff 3 4 145134.0 0.23768 0.00 6.99 8.70 #Na-Cl +pair_coeff 4 4 405774.0 0.23768 0.00 72.40 145.40 #Cl-Cl + +bond_style harmonic +bond_coeff 1 63.014 0.0 +bond_coeff 2 25.724 0.0 + +# ------------------------ Equilibration Run ------------------------------- + +reset_timestep 0 + +thermo 50 +thermo_style custom step etotal pe ke temp press epair evdwl ecoul elong ebond fnorm fmax vol + +compute CSequ all temp/cs cores shells + +# output via chunk method + +#compute prop all property/atom i_CSID +#compute cs_chunk all chunk/atom c_prop +#compute cstherm all temp/chunk cs_chunk temp internal com yes cdof 3.0 +#fix ave_chunk all ave/time 100 1 100 c_cstherm file chunk.dump mode vector + +thermo_modify temp CSequ + +# velocity bias option + +velocity all create 1427 134 dist gaussian mom yes rot no bias yes temp CSequ +Neighbor list info ... + 1 neighbor list requests + update every 1 steps, delay 10 steps, check yes + max neighbors/atom: 2000, page size: 100000 + master list distance cutoff = 22 + ghost atom cutoff = 22 + binsize = 11 -> bins = 3 3 3 +velocity all scale 1427 temp CSequ + +fix thermoberendsen all temp/berendsen 1427 1427 0.4 +fix nve all nve +fix_modify thermoberendsen temp CSequ + +# 2 fmsec timestep + +timestep 0.002 +run 500 +Memory usage per processor = 6.59457 Mbytes +Step TotEng PotEng KinEng Temp Press E_pair E_vdwl E_coul E_long E_bond Fnorm Fmax Volume + 0 -635.80596 -675.46362 39.657659 1427 -21302.622 -675.46362 1.6320365 -677.09565 0 0 1.4528722e-14 2.747802e-15 13990.5 + 50 -634.54091 -666.58532 32.044415 1153.0529 -4531.3775 -668.37478 37.774715 -706.14949 0 1.7894556 9.6796061 2.2085323 13990.5 + 100 -632.30603 -662.37331 30.067272 1081.9095 -3494.2387 -664.52701 39.193627 -703.72064 0 2.1537086 11.074727 2.1546198 13990.5 + 150 -630.53361 -663.40981 32.876196 1182.9829 -73.383109 -665.83886 46.267578 -712.10644 0 2.4290532 11.741311 2.7510336 13990.5 + 200 -628.96388 -664.35527 35.391398 1273.4873 -1706.4631 -666.68771 41.799509 -708.48722 0 2.3324396 10.596222 3.0017211 13990.5 + 250 -627.71591 -662.45411 34.738199 1249.9833 -1256.9123 -665.05391 43.020328 -708.07424 0 2.5998001 10.583081 1.8375441 13990.5 + 300 -627.00338 -664.01808 37.014699 1331.8985 -1884.1512 -666.29225 40.78768 -707.07994 0 2.2741714 9.4246938 1.2451114 13990.5 + 350 -626.21003 -664.96799 38.757968 1394.6265 -1431.9753 -667.39487 41.866515 -709.26138 0 2.426873 10.31352 1.9959989 13990.5 + 400 -625.51909 -661.88115 36.362062 1308.4147 -366.8119 -664.48329 44.342045 -708.82533 0 2.60214 11.071369 2.1790546 13990.5 + 450 -624.9333 -661.08821 36.154909 1300.9607 -2315.3532 -663.62656 39.690513 -703.31708 0 2.5383564 10.013961 1.4015893 13990.5 + 500 -624.02542 -660.36457 36.339145 1307.59 374.84327 -663.00031 45.787267 -708.78758 0 2.6357437 10.67964 2.2503424 13990.5 +Loop time of 5.89773 on 2 procs for 500 steps with 432 atoms + +Performance: 14.650 ns/day, 1.638 hours/ns, 84.778 timesteps/s +99.9% CPU use with 2 MPI tasks x no OpenMP threads + +MPI task timing breakdown: +Section | min time | avg time | max time |%varavg| %total +--------------------------------------------------------------- +Pair | 5.2578 | 5.4112 | 5.5647 | 6.6 | 91.75 +Bond | 0.0015583 | 0.0015992 | 0.0016401 | 0.1 | 0.03 +Neigh | 0.076808 | 0.076815 | 0.076823 | 0.0 | 1.30 +Comm | 0.23842 | 0.39185 | 0.54528 | 24.5 | 6.64 +Output | 0.00029778 | 0.00038707 | 0.00047636 | 0.5 | 0.01 +Modify | 0.0082664 | 0.0083008 | 0.0083351 | 0.0 | 0.14 +Other | | 0.007544 | | | 0.13 + +Nlocal: 216 ave 220 max 212 min +Histogram: 1 0 0 0 0 0 0 0 0 1 +Nghost: 7852 ave 7890 max 7814 min +Histogram: 1 0 0 0 0 0 0 0 0 1 +Neighs: 148601 ave 149750 max 147452 min +Histogram: 1 0 0 0 0 0 0 0 0 1 + +Total # of neighbors = 297202 +Ave neighs/atom = 687.968 +Ave special neighs/atom = 1 +Neighbor list builds = 21 +Dangerous builds = 0 + +unfix thermoberendsen + +# ------------------------ Dynamic Run ------------------------------- + +run 1000 +Memory usage per processor = 6.59743 Mbytes +Step TotEng PotEng KinEng Temp Press E_pair E_vdwl E_coul E_long E_bond Fnorm Fmax Volume + 500 -624.02542 -660.36457 36.339145 1307.59 374.84327 -663.00031 45.787267 -708.78758 0 2.6357437 10.67964 2.2503424 13990.5 + 550 -624.01275 -659.6763 35.663547 1283.28 -408.77541 -662.43528 44.140494 -706.57577 0 2.7589824 11.231423 1.8578347 13990.5 + 600 -624.01608 -661.00741 36.991333 1331.0577 97.833921 -663.55203 45.177756 -708.72979 0 2.5446206 11.407616 1.9340498 13990.5 + 650 -624.01489 -661.39357 37.378678 1344.9955 -1635.0963 -663.98652 41.03269 -705.01921 0 2.5929525 9.7182838 1.8437992 13990.5 + 700 -624.01641 -660.60919 36.592785 1316.7168 274.84427 -663.34199 45.758937 -709.10093 0 2.7327992 11.40835 2.6926812 13990.5 + 750 -624.01423 -660.30947 36.295243 1306.0103 -424.48119 -662.95011 43.863553 -706.81366 0 2.6406337 11.369608 1.7919396 13990.5 + 800 -624.01561 -662.42457 38.408956 1382.068 -1601.0688 -664.8166 41.228621 -706.04522 0 2.3920376 10.379512 1.4876332 13990.5 + 850 -624.01051 -660.5341 36.523587 1314.2268 -1834.5447 -662.93225 40.927147 -703.8594 0 2.3981558 10.724301 2.0240776 13990.5 + 900 -624.01424 -662.4875 38.473259 1384.3818 -285.09099 -665.18917 44.16937 -709.35854 0 2.7016741 11.202878 2.9137876 13990.5 + 950 -624.01787 -661.88077 37.862897 1362.4192 294.61522 -664.50863 45.54898 -710.05761 0 2.6278644 12.347191 2.5194653 13990.5 + 1000 -624.01006 -658.82908 34.819022 1252.8915 1297.8621 -661.39405 48.235238 -709.62929 0 2.5649692 13.763173 2.6809836 13990.5 + 1050 -624.01237 -657.96252 33.950153 1221.627 -565.05956 -661.022 43.831206 -704.85321 0 3.059482 10.980185 1.5069363 13990.5 + 1100 -624.01221 -662.32641 38.314199 1378.6583 -1053.4531 -664.68928 42.37808 -707.06736 0 2.3628693 11.058878 2.310501 13990.5 + 1150 -624.01355 -663.31571 39.30216 1414.2081 -682.32516 -665.67212 42.824978 -708.4971 0 2.3564095 9.9617506 1.4794302 13990.5 + 1200 -624.01228 -660.54648 36.534208 1314.609 -1922.8104 -663.09623 40.35503 -703.45126 0 2.5497509 10.105408 1.5577027 13990.5 + 1250 -624.01453 -660.19988 36.18535 1302.056 1030.8826 -662.54352 47.329262 -709.87278 0 2.3436398 11.480691 2.3100697 13990.5 + 1300 -624.0174 -660.44812 36.430711 1310.8848 -592.37472 -663.42007 43.498393 -706.91846 0 2.9719516 10.517346 1.4566735 13990.5 + 1350 -624.02147 -662.7255 38.704027 1392.6855 328.82353 -665.26912 45.2889 -710.55802 0 2.5436194 11.071493 2.299515 13990.5 + 1400 -624.01155 -660.32249 36.310942 1306.5752 664.45071 -663.12835 46.391793 -709.52014 0 2.8058565 12.0411 2.0813357 13990.5 + 1450 -624.01128 -660.37598 36.364706 1308.5098 -954.57186 -662.9841 42.653305 -705.63741 0 2.6081222 10.498899 1.7930019 13990.5 + 1500 -624.01495 -658.96253 34.947585 1257.5176 14.253508 -661.8179 45.019135 -706.83704 0 2.8553699 10.911364 1.8897264 13990.5 +Loop time of 11.7651 on 2 procs for 1000 steps with 432 atoms + +Performance: 14.687 ns/day, 1.634 hours/ns, 84.997 timesteps/s +99.9% CPU use with 2 MPI tasks x no OpenMP threads + +MPI task timing breakdown: +Section | min time | avg time | max time |%varavg| %total +--------------------------------------------------------------- +Pair | 10.824 | 10.886 | 10.948 | 1.9 | 92.53 +Bond | 0.0032091 | 0.0032401 | 0.0032711 | 0.1 | 0.03 +Neigh | 0.16206 | 0.16206 | 0.16206 | 0.0 | 1.38 +Comm | 0.62978 | 0.6916 | 0.75341 | 7.4 | 5.88 +Output | 0.00071597 | 0.00087166 | 0.0010273 | 0.5 | 0.01 +Modify | 0.0043559 | 0.0043788 | 0.0044017 | 0.0 | 0.04 +Other | | 0.01691 | | | 0.14 + +Nlocal: 216 ave 224 max 208 min +Histogram: 1 0 0 0 0 0 0 0 0 1 +Nghost: 7910 ave 7994 max 7826 min +Histogram: 1 0 0 0 0 0 0 0 0 1 +Neighs: 148484 ave 151739 max 145228 min +Histogram: 1 0 0 0 0 0 0 0 0 1 + +Total # of neighbors = 296967 +Ave neighs/atom = 687.424 +Ave special neighs/atom = 1 +Neighbor list builds = 46 +Dangerous builds = 0 +Total wall time: 0:00:17 diff --git a/src/.gitignore b/src/.gitignore index 632b46a08a..02b05b364b 100644 --- a/src/.gitignore +++ b/src/.gitignore @@ -976,6 +976,8 @@ /fix_ttm_mod.h /pair_born_coul_long_cs.cpp /pair_born_coul_long_cs.h +/pair_born_coul_dsf_cs.cpp +/pair_born_coul_dsf_cs.h /pair_buck_coul_long_cs.cpp /pair_buck_coul_long_cs.h /pair_coul_long_cs.cpp diff --git a/src/CORESHELL/Install.sh b/src/CORESHELL/Install.sh index f5ea54ac86..7c0b7a02a2 100644 --- a/src/CORESHELL/Install.sh +++ b/src/CORESHELL/Install.sh @@ -31,8 +31,10 @@ action () { action compute_temp_cs.cpp action compute_temp_cs.h action pair_born_coul_long_cs.cpp pair_born_coul_long.cpp +action pair_born_coul_dsf_cs.cpp pair_born_coul_dsf.cpp action pair_buck_coul_long_cs.cpp pair_buck_coul_long.cpp action pair_born_coul_long_cs.h pair_born_coul_long.h +action pair_born_coul_dsf_cs.h pair_born_coul_dsf.h action pair_buck_coul_long_cs.h pair_buck_coul_long.h action pair_coul_long_cs.cpp pair_coul_long.cpp action pair_coul_long_cs.h pair_coul_long.h diff --git a/src/CORESHELL/pair_born_coul_dsf_cs.cpp b/src/CORESHELL/pair_born_coul_dsf_cs.cpp new file mode 100644 index 0000000000..2cd2f211ff --- /dev/null +++ b/src/CORESHELL/pair_born_coul_dsf_cs.cpp @@ -0,0 +1,162 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Ariel Lozano (arielzn@gmail.com) + References: Fennell and Gezelter, JCP 124, 234104 (2006) +------------------------------------------------------------------------- */ + +#include +#include +#include +#include +#include "pair_born_coul_dsf_cs.h" +#include "atom.h" +#include "comm.h" +#include "force.h" +#include "neighbor.h" +#include "neigh_list.h" +#include "math_const.h" +#include "memory.h" +#include "error.h" +#include "math_special.h" + +using namespace LAMMPS_NS; +using namespace MathConst; + +#define EPSILON 1.0e-20 + + +/* ---------------------------------------------------------------------- */ + +PairBornCoulDSFCS::PairBornCoulDSFCS(LAMMPS *lmp) : PairBornCoulDSF(lmp) +{ + writedata = 1; + single_enable = 0; +} + +/* ---------------------------------------------------------------------- */ + +void PairBornCoulDSFCS::compute(int eflag, int vflag) +{ + int i,j,ii,jj,inum,jnum,itype,jtype; + double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair; + double r,rsq,r2inv,r6inv,forcecoul,forceborn,factor_coul,factor_lj; + double prefactor,erfcc,erfcd,arg; + double rexp; + int *ilist,*jlist,*numneigh,**firstneigh; + + evdwl = ecoul = 0.0; + if (eflag || vflag) ev_setup(eflag,vflag); + else evflag = vflag_fdotr = 0; + + double **x = atom->x; + double **f = atom->f; + double *q = atom->q; + int *type = atom->type; + int nlocal = atom->nlocal; + double *special_lj = force->special_lj; + double *special_coul = force->special_coul; + int newton_pair = force->newton_pair; + double qqrd2e = force->qqrd2e; + + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + + // loop over neighbors of my atoms + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + qtmp = q[i]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + itype = type[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + + // self coulombic energy + if (eflag) { + double e_self = -(e_shift/2.0 + alpha/MY_PIS) * qtmp*qtmp*qqrd2e; + ev_tally(i,i,nlocal,0,0.0,e_self,0.0,0.0,0.0,0.0); + } + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + factor_lj = special_lj[sbmask(j)]; + factor_coul = special_coul[sbmask(j)]; + j &= NEIGHMASK; + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx*delx + dely*dely + delz*delz; + jtype = type[j]; + + if (rsq < cutsq[itype][jtype]) { + rsq += EPSILON; // Add Epsilon for case: r = 0; Interaction must be removed by special bond; + r2inv = 1.0/rsq; + + if (rsq < cut_coulsq) { + r = sqrt(rsq); + prefactor = qqrd2e*qtmp*q[j] / r; + arg = alpha * r ; + erfcd = MathSpecial::expmsq(arg); + erfcc = MathSpecial::my_erfcx(arg) * erfcd; + forcecoul = prefactor * (erfcc/r + 2.0*alpha/MY_PIS * erfcd + + r*f_shift) * r; + if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor; + } else forcecoul = 0.0; + + if (rsq < cut_ljsq[itype][jtype]) { + r6inv = r2inv*r2inv*r2inv; + r = sqrt(rsq); + rexp = exp((sigma[itype][jtype]-r)*rhoinv[itype][jtype]); + forceborn = born1[itype][jtype]*r*rexp - born2[itype][jtype]*r6inv + + born3[itype][jtype]*r2inv*r6inv; + } else forceborn = 0.0; + + fpair = (forcecoul + factor_lj*forceborn) * r2inv; + + f[i][0] += delx*fpair; + f[i][1] += dely*fpair; + f[i][2] += delz*fpair; + if (newton_pair || j < nlocal) { + f[j][0] -= delx*fpair; + f[j][1] -= dely*fpair; + f[j][2] -= delz*fpair; + } + + if (eflag) { + if (rsq < cut_coulsq) { + ecoul = prefactor * (erfcc - r*e_shift - rsq*f_shift); + if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor; + } else ecoul = 0.0; + if (rsq < cut_ljsq[itype][jtype]) { + evdwl = a[itype][jtype]*rexp - c[itype][jtype]*r6inv + + d[itype][jtype]*r6inv*r2inv - offset[itype][jtype]; + evdwl *= factor_lj; + } else evdwl = 0.0; + } + + if (evflag) ev_tally(i,j,nlocal,newton_pair, + evdwl,ecoul,fpair,delx,dely,delz); + } + } + } + + if (vflag_fdotr) virial_fdotr_compute(); +} + diff --git a/src/CORESHELL/pair_born_coul_dsf_cs.h b/src/CORESHELL/pair_born_coul_dsf_cs.h new file mode 100644 index 0000000000..a49df5971c --- /dev/null +++ b/src/CORESHELL/pair_born_coul_dsf_cs.h @@ -0,0 +1,59 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef PAIR_CLASS + +PairStyle(born/coul/dsf/cs,PairBornCoulDSFCS) + +#else + +#ifndef LMP_PAIR_BORN_COUL_DSF_CS_H +#define LMP_PAIR_BORN_COUL_DSF_CS_H + +#include "pair_born_coul_dsf.h" + +namespace LAMMPS_NS { + +class PairBornCoulDSFCS : public PairBornCoulDSF { + public: + PairBornCoulDSFCS(class LAMMPS *); + virtual void compute(int, int); +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +E: Illegal ... command + +Self-explanatory. Check the input script syntax and compare to the +documentation for the command. You can use -echo screen as a +command-line option when running LAMMPS to see the offending line. + +E: Incorrect args for pair coefficients + +Self-explanatory. Check the input script or data file. + +E: Pair style born/coul/dsf requires atom attribute q + +The atom style defined does not have this attribute. + +E: All pair coeffs are not set + +All pair coefficients must be set in the data file or by the +pair_coeff command before running a simulation. + +*/ diff --git a/src/pair_born_coul_dsf.cpp b/src/pair_born_coul_dsf.cpp new file mode 100644 index 0000000000..a79790dc48 --- /dev/null +++ b/src/pair_born_coul_dsf.cpp @@ -0,0 +1,494 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Ariel Lozano (arielzn@gmail.com) + References: Fennell and Gezelter, JCP 124, 234104 (2006) +------------------------------------------------------------------------- */ + +#include +#include +#include +#include +#include "pair_born_coul_dsf.h" +#include "atom.h" +#include "comm.h" +#include "force.h" +#include "neighbor.h" +#include "neigh_list.h" +#include "math_const.h" +#include "memory.h" +#include "error.h" +#include "math_special.h" + +using namespace LAMMPS_NS; +using namespace MathConst; + +/* ---------------------------------------------------------------------- */ + +PairBornCoulDSF::PairBornCoulDSF(LAMMPS *lmp) : Pair(lmp) +{ + writedata = 1; + single_enable = 0; +} + +/* ---------------------------------------------------------------------- */ + +PairBornCoulDSF::~PairBornCoulDSF() +{ + if (allocated) { + memory->destroy(setflag); + memory->destroy(cutsq); + + memory->destroy(cut_lj); + memory->destroy(cut_ljsq); + memory->destroy(a); + memory->destroy(rho); + memory->destroy(sigma); + memory->destroy(c); + memory->destroy(d); + memory->destroy(rhoinv); + memory->destroy(born1); + memory->destroy(born2); + memory->destroy(born3); + memory->destroy(offset); + } +} + +/* ---------------------------------------------------------------------- */ + +void PairBornCoulDSF::compute(int eflag, int vflag) +{ + int i,j,ii,jj,inum,jnum,itype,jtype; + double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair; + double r,rsq,r2inv,r6inv,forcecoul,forceborn,factor_coul,factor_lj; + double prefactor,erfcc,erfcd,arg; + double rexp; + int *ilist,*jlist,*numneigh,**firstneigh; + + evdwl = ecoul = 0.0; + if (eflag || vflag) ev_setup(eflag,vflag); + else evflag = vflag_fdotr = 0; + + double **x = atom->x; + double **f = atom->f; + double *q = atom->q; + int *type = atom->type; + int nlocal = atom->nlocal; + double *special_lj = force->special_lj; + double *special_coul = force->special_coul; + int newton_pair = force->newton_pair; + double qqrd2e = force->qqrd2e; + + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + + // loop over neighbors of my atoms + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + qtmp = q[i]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + itype = type[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + + // self coulombic energy + if (eflag) { + double e_self = -(e_shift/2.0 + alpha/MY_PIS) * qtmp*qtmp*qqrd2e; + ev_tally(i,i,nlocal,0,0.0,e_self,0.0,0.0,0.0,0.0); + } + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + factor_lj = special_lj[sbmask(j)]; + factor_coul = special_coul[sbmask(j)]; + j &= NEIGHMASK; + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx*delx + dely*dely + delz*delz; + jtype = type[j]; + + if (rsq < cutsq[itype][jtype]) { + r2inv = 1.0/rsq; + + if (rsq < cut_coulsq) { + r = sqrt(rsq); + prefactor = qqrd2e*qtmp*q[j]/r; + arg = alpha * r ; + erfcd = MathSpecial::expmsq(arg); + erfcc = MathSpecial::my_erfcx(arg) * erfcd; + forcecoul = prefactor * (erfcc/r + 2.0*alpha/MY_PIS * erfcd + + r*f_shift) * r; + if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor; + } else forcecoul = 0.0; + + if (rsq < cut_ljsq[itype][jtype]) { + r6inv = r2inv*r2inv*r2inv; + r = sqrt(rsq); + rexp = exp((sigma[itype][jtype]-r)*rhoinv[itype][jtype]); + forceborn = born1[itype][jtype]*r*rexp - born2[itype][jtype]*r6inv + + born3[itype][jtype]*r2inv*r6inv; + } else forceborn = 0.0; + + fpair = (forcecoul + factor_lj*forceborn) * r2inv; + + f[i][0] += delx*fpair; + f[i][1] += dely*fpair; + f[i][2] += delz*fpair; + if (newton_pair || j < nlocal) { + f[j][0] -= delx*fpair; + f[j][1] -= dely*fpair; + f[j][2] -= delz*fpair; + } + + if (eflag) { + if (rsq < cut_coulsq) { + ecoul = prefactor * (erfcc - r*e_shift - rsq*f_shift); + if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor; + } else ecoul = 0.0; + if (rsq < cut_ljsq[itype][jtype]) { + evdwl = a[itype][jtype]*rexp - c[itype][jtype]*r6inv + + d[itype][jtype]*r6inv*r2inv - offset[itype][jtype]; + evdwl *= factor_lj; + } else evdwl = 0.0; + } + + if (evflag) ev_tally(i,j,nlocal,newton_pair, + evdwl,ecoul,fpair,delx,dely,delz); + } + } + } + + if (vflag_fdotr) virial_fdotr_compute(); +} + +/* ---------------------------------------------------------------------- + allocate all arrays +------------------------------------------------------------------------- */ + +void PairBornCoulDSF::allocate() +{ + allocated = 1; + int n = atom->ntypes; + + memory->create(setflag,n+1,n+1,"pair:setflag"); + for (int i = 1; i <= n; i++) + for (int j = i; j <= n; j++) + setflag[i][j] = 0; + + memory->create(cutsq,n+1,n+1,"pair:cutsq"); + + memory->create(cut_lj,n+1,n+1,"pair:cut_lj"); + memory->create(cut_ljsq,n+1,n+1,"pair:cut_ljsq"); + memory->create(a,n+1,n+1,"pair:a"); + memory->create(rho,n+1,n+1,"pair:rho"); + memory->create(sigma,n+1,n+1,"pair:sigma"); + memory->create(c,n+1,n+1,"pair:c"); + memory->create(d,n+1,n+1,"pair:d"); + memory->create(rhoinv,n+1,n+1,"pair:rhoinv"); + memory->create(born1,n+1,n+1,"pair:born1"); + memory->create(born2,n+1,n+1,"pair:born2"); + memory->create(born3,n+1,n+1,"pair:born3"); + memory->create(offset,n+1,n+1,"pair:offset"); +} + +/* ---------------------------------------------------------------------- + global settings +------------------------------------------------------------------------- */ + +void PairBornCoulDSF::settings(int narg, char **arg) +{ + if (narg < 2 || narg > 3) error->all(FLERR,"Illegal pair_style command"); + + alpha = force->numeric(FLERR,arg[0]); + cut_lj_global = force->numeric(FLERR,arg[1]); + if (narg == 2) cut_coul = cut_lj_global; + else cut_coul = force->numeric(FLERR,arg[2]); + + if (allocated) { + int i,j; + for (i = 1; i <= atom->ntypes; i++) + for (j = i+1; j <= atom->ntypes; j++) + if (setflag[i][j]) + cut_lj[i][j] = cut_lj_global; + } +} + +/* ---------------------------------------------------------------------- + set coeffs for one or more type pairs +------------------------------------------------------------------------- */ + +void PairBornCoulDSF::coeff(int narg, char **arg) +{ + if (narg < 7 || narg > 8) + error->all(FLERR,"Incorrect args for pair coefficients"); + if (!allocated) allocate(); + + int ilo,ihi,jlo,jhi; + force->bounds(FLERR,arg[0],atom->ntypes,ilo,ihi); + force->bounds(FLERR,arg[1],atom->ntypes,jlo,jhi); + + double a_one = force->numeric(FLERR,arg[2]); + double rho_one = force->numeric(FLERR,arg[3]); + double sigma_one = force->numeric(FLERR,arg[4]); + if (rho_one <= 0) error->all(FLERR,"Incorrect args for pair coefficients"); + double c_one = force->numeric(FLERR,arg[5]); + double d_one = force->numeric(FLERR,arg[6]); + + double cut_lj_one = cut_lj_global; + if (narg == 8) cut_lj_one = force->numeric(FLERR,arg[7]); + + int count = 0; + for (int i = ilo; i <= ihi; i++) { + for (int j = MAX(jlo,i); j <= jhi; j++) { + a[i][j] = a_one; + rho[i][j] = rho_one; + sigma[i][j] = sigma_one; + c[i][j] = c_one; + d[i][j] = d_one; + cut_lj[i][j] = cut_lj_one; + setflag[i][j] = 1; + count++; + } + } + + if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); +} + +/* ---------------------------------------------------------------------- + init specific to this pair style +------------------------------------------------------------------------- */ + +void PairBornCoulDSF::init_style() +{ + if (!atom->q_flag) + error->all(FLERR,"Pair style born/coul/dsf requires atom attribute q"); + + neighbor->request(this,instance_me); + + cut_coulsq = cut_coul * cut_coul; + double erfcc = erfc(alpha*cut_coul); + double erfcd = exp(-alpha*alpha*cut_coul*cut_coul); + f_shift = -(erfcc/cut_coulsq + 2.0/MY_PIS*alpha*erfcd/cut_coul); + e_shift = erfcc/cut_coul - f_shift*cut_coul; +} + +/* ---------------------------------------------------------------------- + init for one type pair i,j and corresponding j,i +------------------------------------------------------------------------- */ + +double PairBornCoulDSF::init_one(int i, int j) +{ + if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set"); + + double cut = MAX(cut_lj[i][j],cut_coul); + cut_ljsq[i][j] = cut_lj[i][j] * cut_lj[i][j]; + + rhoinv[i][j] = 1.0/rho[i][j]; + born1[i][j] = a[i][j]/rho[i][j]; + born2[i][j] = 6.0*c[i][j]; + born3[i][j] = 8.0*d[i][j]; + + if (offset_flag) { + double rexp = exp((sigma[i][j]-cut_lj[i][j])*rhoinv[i][j]); + offset[i][j] = a[i][j]*rexp - c[i][j]/pow(cut_lj[i][j],6.0) + + d[i][j]/pow(cut_lj[i][j],8.0); + } else offset[i][j] = 0.0; + + cut_ljsq[j][i] = cut_ljsq[i][j]; + a[j][i] = a[i][j]; + c[j][i] = c[i][j]; + d[j][i] = d[i][j]; + rhoinv[j][i] = rhoinv[i][j]; + sigma[j][i] = sigma[i][j]; + born1[j][i] = born1[i][j]; + born2[j][i] = born2[i][j]; + born3[j][i] = born3[i][j]; + offset[j][i] = offset[i][j]; + + return cut; +} + +/* ---------------------------------------------------------------------- + proc 0 writes to restart file +------------------------------------------------------------------------- */ + +void PairBornCoulDSF::write_restart(FILE *fp) +{ + write_restart_settings(fp); + + int i,j; + for (i = 1; i <= atom->ntypes; i++) + for (j = i; j <= atom->ntypes; j++) { + fwrite(&setflag[i][j],sizeof(int),1,fp); + if (setflag[i][j]) { + fwrite(&a[i][j],sizeof(double),1,fp); + fwrite(&rho[i][j],sizeof(double),1,fp); + fwrite(&sigma[i][j],sizeof(double),1,fp); + fwrite(&c[i][j],sizeof(double),1,fp); + fwrite(&d[i][j],sizeof(double),1,fp); + fwrite(&cut_lj[i][j],sizeof(double),1,fp); + } + } +} + +/* ---------------------------------------------------------------------- + proc 0 reads from restart file, bcasts +------------------------------------------------------------------------- */ + +void PairBornCoulDSF::read_restart(FILE *fp) +{ + read_restart_settings(fp); + + allocate(); + + int i,j; + int me = comm->me; + for (i = 1; i <= atom->ntypes; i++) + for (j = i; j <= atom->ntypes; j++) { + if (me == 0) fread(&setflag[i][j],sizeof(int),1,fp); + MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world); + if (setflag[i][j]) { + if (me == 0) { + fread(&a[i][j],sizeof(double),1,fp); + fread(&rho[i][j],sizeof(double),1,fp); + fread(&sigma[i][j],sizeof(double),1,fp); + fread(&c[i][j],sizeof(double),1,fp); + fread(&d[i][j],sizeof(double),1,fp); + fread(&cut_lj[i][j],sizeof(double),1,fp); + } + MPI_Bcast(&a[i][j],1,MPI_DOUBLE,0,world); + MPI_Bcast(&rho[i][j],1,MPI_DOUBLE,0,world); + MPI_Bcast(&sigma[i][j],1,MPI_DOUBLE,0,world); + MPI_Bcast(&c[i][j],1,MPI_DOUBLE,0,world); + MPI_Bcast(&d[i][j],1,MPI_DOUBLE,0,world); + MPI_Bcast(&cut_lj[i][j],1,MPI_DOUBLE,0,world); + } + } +} + +/* ---------------------------------------------------------------------- + proc 0 writes to restart file +------------------------------------------------------------------------- */ + +void PairBornCoulDSF::write_restart_settings(FILE *fp) +{ + fwrite(&alpha,sizeof(double),1,fp); + fwrite(&cut_lj_global,sizeof(double),1,fp); + fwrite(&cut_coul,sizeof(double),1,fp); + fwrite(&offset_flag,sizeof(int),1,fp); + fwrite(&mix_flag,sizeof(int),1,fp); +} + +/* ---------------------------------------------------------------------- + proc 0 reads from restart file, bcasts +------------------------------------------------------------------------- */ + +void PairBornCoulDSF::read_restart_settings(FILE *fp) +{ + if (comm->me == 0) { + fread(&alpha,sizeof(double),1,fp); + fread(&cut_lj_global,sizeof(double),1,fp); + fread(&cut_coul,sizeof(double),1,fp); + fread(&offset_flag,sizeof(int),1,fp); + fread(&mix_flag,sizeof(int),1,fp); + } + MPI_Bcast(&alpha,1,MPI_DOUBLE,0,world); + MPI_Bcast(&cut_lj_global,1,MPI_DOUBLE,0,world); + MPI_Bcast(&cut_coul,1,MPI_DOUBLE,0,world); + MPI_Bcast(&offset_flag,1,MPI_INT,0,world); + MPI_Bcast(&mix_flag,1,MPI_INT,0,world); +} + +/* ---------------------------------------------------------------------- + proc 0 writes to data file +------------------------------------------------------------------------- */ + +void PairBornCoulDSF::write_data(FILE *fp) +{ + for (int i = 1; i <= atom->ntypes; i++) + fprintf(fp,"%d %g %g %g %g %g\n",i, + a[i][i],rho[i][i],sigma[i][i],c[i][i],d[i][i]); +} + +/* ---------------------------------------------------------------------- + proc 0 writes all pairs to data file +------------------------------------------------------------------------- */ + +void PairBornCoulDSF::write_data_all(FILE *fp) +{ + for (int i = 1; i <= atom->ntypes; i++) + for (int j = i; j <= atom->ntypes; j++) + fprintf(fp,"%d %d %g %g %g %g %g %g\n",i,j, + a[i][j],rho[i][j],sigma[i][j],c[i][j],d[i][j],cut_lj[i][j]); +} + +/* ---------------------------------------------------------------------- + only the pair part is calculated here +------------------------------------------------------------------------- */ + +double PairBornCoulDSF::single(int i, int j, int itype, int jtype, + double rsq, + double factor_coul, double factor_lj, + double &fforce) +{ + double r2inv,r6inv,r,prefactor,rexp; + double forcecoul,forceborn,phicoul,phiborn; + double erfcc,erfcd,arg; + + r2inv = 1.0/rsq; + + if (rsq < cut_coulsq) { + r = sqrt(rsq); + prefactor = factor_coul * force->qqrd2e * atom->q[i]*atom->q[j]/r; + + arg = alpha * r ; + erfcd = MathSpecial::expmsq(arg); + erfcc = MathSpecial::my_erfcx(arg) * erfcd; + + forcecoul = prefactor * (erfcc/r + 2.0*alpha/MY_PIS*erfcd + + r*f_shift) * r; + + if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor; + } else forcecoul = 0.0; + + if (rsq < cut_ljsq[itype][jtype]) { + r6inv = r2inv*r2inv*r2inv; + r = sqrt(rsq); + rexp = exp(-r*rhoinv[itype][jtype]); + forceborn = born1[itype][jtype]*r*rexp - born2[itype][jtype]*r6inv + + born3[itype][jtype]*r2inv*r6inv; + } else forceborn = 0.0; + + fforce = (forcecoul + factor_lj*forceborn) * r2inv; + + double eng = 0.0; + if (rsq < cut_coulsq) { + phicoul = prefactor * (erfcc - r*e_shift - rsq*f_shift); + if (factor_coul < 1.0) phicoul -= (1.0-factor_coul)*prefactor; + eng += phicoul; + } + if (rsq < cut_ljsq[itype][jtype]) { + phiborn = a[itype][jtype]*rexp - c[itype][jtype]*r6inv + + d[itype][jtype]*r2inv*r6inv - offset[itype][jtype]; + eng += factor_lj*phiborn; + } + return eng; +} diff --git a/src/pair_born_coul_dsf.h b/src/pair_born_coul_dsf.h new file mode 100644 index 0000000000..c7cfb26094 --- /dev/null +++ b/src/pair_born_coul_dsf.h @@ -0,0 +1,81 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef PAIR_CLASS + +PairStyle(born/coul/dsf,PairBornCoulDSF) + +#else + +#ifndef LMP_PAIR_BORN_COUL_DSF_H +#define LMP_PAIR_BORN_COUL_DSF_H + +#include "pair.h" + +namespace LAMMPS_NS { + +class PairBornCoulDSF : public Pair { + public: + PairBornCoulDSF(class LAMMPS *); + virtual ~PairBornCoulDSF(); + virtual void compute(int, int); + void settings(int, char **); + void coeff(int, char **); + void init_style(); + double init_one(int, int); + void write_restart(FILE *); + void read_restart(FILE *); + void write_restart_settings(FILE *); + void read_restart_settings(FILE *); + void write_data(FILE *); + void write_data_all(FILE *); + double single(int, int, int, int, double, double, double, double &); + + protected: + double cut_lj_global,alpha; + double **cut_lj,**cut_ljsq; + double cut_coul,cut_coulsq; + double **a,**rho,**sigma,**c,**d; + double **rhoinv,**born1,**born2,**born3,**offset; + double f_shift,e_shift; + + void allocate(); +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +E: Illegal ... command + +Self-explanatory. Check the input script syntax and compare to the +documentation for the command. You can use -echo screen as a +command-line option when running LAMMPS to see the offending line. + +E: Incorrect args for pair coefficients + +Self-explanatory. Check the input script or data file. + +E: Pair style born/coul/dsf requires atom attribute q + +The atom style defined does not have this attribute. + +E: All pair coeffs are not set + +All pair coefficients must be set in the data file or by the +pair_coeff command before running a simulation. + +*/