diff --git a/doc/src/fix_eos_table_rx.txt b/doc/src/fix_eos_table_rx.txt
index f92b405f495b07a9bf690f8aa0aaafeb23f65813..749642f57c203c753bc51fc425be0b56fed5b6a7 100644
--- a/doc/src/fix_eos_table_rx.txt
+++ b/doc/src/fix_eos_table_rx.txt
@@ -10,7 +10,7 @@ fix eos/table/rx command :h3
[Syntax:]
-fix ID group-ID eos/table/rx style file1 N keyword file2 :pre
+fix ID group-ID eos/table/rx style file1 N keyword ... :pre
ID, group-ID are documented in "fix"_fix.html command
eos/table/rx = style name of this fix command
@@ -18,11 +18,16 @@ style = {linear} = method of interpolation
file1 = filename containing the tabulated equation of state
N = use N values in {linear} tables
keyword = name of table keyword correponding to table file
-file2 = filename containing the heats of formation of each species :ul
+file2 = filename containing the heats of formation of each species (optional)
+deltaHf = heat of formation for a single species in energy units (optional)
+energyCorr = energy correction in energy units (optional)
+tempCorrCoeff = temperature correction coefficient (optional) :ul
[Examples:]
-fix 1 all eos/table/rx linear eos.table 10000 KEYWORD thermo.table :pre
+fix 1 all eos/table/rx linear eos.table 10000 KEYWORD thermo.table
+fix 1 all eos/table/rx linear eos.table 10000 KEYWORD 1.5
+fix 1 all eos/table/rx linear eos.table 10000 KEYWORD 1.5 0.025 0.0 :pre
[Description:]
@@ -39,7 +44,15 @@ where {m} is the number of species, {c_i,j} is the concentration of
species {j} in particle {i}, {u_j} is the internal energy of species j,
{DeltaH_f,j} is the heat of formation of species {j}, N is the number of
molecules represented by the coarse-grained particle, kb is the
-Boltzmann constant, and T is the temperature of the system.
+Boltzmann constant, and T is the temperature of the system. Additionally,
+it is possible to modify the concentration-dependent particle internal
+energy relation by adding an energy correction, temperature-dependent
+correction, and/or a molecule-dependent correction. An energy correction can
+be specified as a constant (in energy units). A temperature correction can be
+specified by multiplying a temperature correction coefficient by the
+internal temperature. A molecular correction can be specified by
+by multiplying a molecule correction coefficient by the average number of
+product gas particles in the coarse-grain particle.
Fix {eos/table/rx} creates interpolation tables of length {N} from {m}
internal energy values of each species {u_j} listed in a file as a
@@ -58,6 +71,14 @@ file is described below.
The second filename specifies a file containing heat of formation
{DeltaH_f,j} for each species.
+In cases where the coarse-grain particle represents a single molecular
+species (i.e., no reactions occur and fix {rx} is not present in the input file),
+fix {eos/table/rx} can be applied in a similar manner to fix {eos/table}
+within a non-reactive DPD simulation. In this case, the heat of formation
+filename is replaced with the heat of formation value for the single species.
+Additionally, the energy correction and temperature correction coefficients may
+also be specified as fix arguments.
+
:line
The format of a tabulated file is as follows (without the
@@ -116,6 +137,19 @@ Note that the species can be listed in any order. The tag that is
used as the species name must correspond with the tags used to define
the reactions with the "fix rx"_fix_rx.html command.
+Alternatively, corrections to the EOS can be included by specifying
+three additional columns that correspond to the energy correction,
+the temperature correction coefficient and molecule correction
+coefficient. In this case, the format of the file is as follows:
+
+# HEAT OF FORMATION TABLE (one or more comment or blank lines) :pre
+ (blank)
+h2 0.00 1.23 0.025 0.0 (species name, heat of formation, energy correction, temperature correction coefficient, molecule correction coefficient)
+no2 0.34 0.00 0.000 -1.76
+n2 0.00 0.00 0.000 -1.76
+...
+no 0.93 0.00 0.000 -1.76 :pre
+
:line
[Restrictions:]
diff --git a/doc/src/pair_exp6_rx.txt b/doc/src/pair_exp6_rx.txt
index 7b22dccc4fb49448d9136f0df3411ca0c2957ac3..dafba2c44c09f34e92b66f0edcba6db7889a8f21 100644
--- a/doc/src/pair_exp6_rx.txt
+++ b/doc/src/pair_exp6_rx.txt
@@ -10,16 +10,21 @@ pair_style exp6/rx command :h3
[Syntax:]
-pair_style exp6/rx cutoff :pre
+pair_style exp6/rx cutoff ... :pre
-cutoff = global cutoff for DPD interactions (distance units) :ul
+cutoff = global cutoff for DPD interactions (distance units)
+weighting = fractional or molecular (optional) :ul
[Examples:]
pair_style exp6/rx 10.0
-pair_coeff * * exp6.params h2o h2o 1.0 1.0 10.0
-pair_coeff * * exp6.params h2o 1fluid 1.0 1.0 10.0
-pair_coeff * * exp6.params 1fluid 1fluid 1.0 1.0 10.0 :pre
+pair_style exp6/rx 10.0 fractional
+pair_style exp6/rx 10.0 molecular
+pair_coeff * * exp6.params h2o h2o exponent 1.0 1.0 10.0
+pair_coeff * * exp6.params h2o 1fluid exponent 1.0 1.0 10.0
+pair_coeff * * exp6.params 1fluid 1fluid exponent 1.0 1.0 10.0
+pair_coeff * * exp6.params 1fluid 1fluid none 10.0
+pair_coeff * * exp6.params 1fluid 1fluid polynomial filename 10.0 :pre
[Description:]
@@ -50,14 +55,36 @@ defined in the reaction kinetics files specified with the "fix
rx"_fix_rx.html command or they must correspond to the tag "1fluid",
signifying interaction with a product species mixture determined
through a one-fluid approximation. The interaction potential is
-weighted by the geometric average of the concentrations of the two
-species. The coarse-grained potential is stored before and after the
+weighted by the geometric average of either the mole fraction concentrations
+or the number of molecules associated with the interacting coarse-grained
+particles (see the {fractional} or {molecular} weighting pair style options).
+The coarse-grained potential is stored before and after the
reaction kinetics solver is applied, where the difference is defined
to be the internal chemical energy (uChem).
-The fourth and fifth arguments specify the {Rm} and {epsilon} scaling exponents.
+The fourth argument specifies the type of scaling that will be used
+to scale the EXP-6 paramters as reactions occur. Currently, there
+are three scaling options: {exponent}, {polynomial} and {none}.
-The final argument specifies the interaction cutoff.
+Exponent scaling requires two additional arguments for scaling
+the {Rm} and {epsilon} parameters, respectively. The scaling factor
+is computed by phi^exponent, where phi is the number of molecules
+represented by the coarse-grain particle and exponent is specified
+as a pair coefficient argument for {Rm} and {epsilon}, respectively.
+The {Rm} and {epsilon} parameters are multiplied by the scaling
+factor to give the scaled interaction paramters for the CG particle.
+
+Polynomial scaling requires a filename to be specified as a pair
+coeff argument. The file contains the coefficients to a fifth order
+polynomial for the {alpha}, {epsilon} and {Rm} parameters that depend
+upon phi (the number of molecules represented by the CG particle).
+The format of a polynomial file is provided below.
+
+The {none} option to the scaling does not have any additional pair coeff
+arguments. This is equivalent to specifying the {exponent} option with
+{Rm} and {epsilon} exponents of 0.0 and 0.0, respectively.
+
+The final argument specifies the interaction cutoff (optional).
:line
@@ -70,6 +97,19 @@ no2 exp6 13.60 0.01 3.70
...
co2 exp6 13.00 0.03 3.20 :pre
+The format of the polynomial scaling file as follows (without the
+parenthesized comments):
+
+# POLYNOMIAL FILE (one or more comment or blank lines) :pre
+# General Functional Form:
+# A*phi^5 + B*phi^4 + C*phi^3 + D*phi^2 + E*phi + F
+#
+# Parameter A B C D E F
+ (blank)
+alpha 0.0000 0.00000 0.00008 0.04955 -0.73804 13.63201
+epsilon 0.0000 0.00478 -0.06283 0.24486 -0.33737 2.60097
+rm 0.0001 -0.00118 -0.00253 0.05812 -0.00509 1.50106 :pre
+
A section begins with a non-blank line whose 1st character is not a
"#"; blank lines or lines starting with "#" can be used as comments
between sections.
@@ -117,4 +157,4 @@ LAMMPS"_Section_start.html#start_3 section for more info.
"pair_coeff"_pair_coeff.html
-[Default:] none
+[Default:] fractional weighting
diff --git a/doc/src/pair_multi_lucy_rx.txt b/doc/src/pair_multi_lucy_rx.txt
index 14b5b32181e6b561a0bf8f602b446e17c0b32eb7..75547a71cedef8c3d910d722c2edc90188b6d5f1 100644
--- a/doc/src/pair_multi_lucy_rx.txt
+++ b/doc/src/pair_multi_lucy_rx.txt
@@ -13,11 +13,14 @@ pair_style multi/lucy/rx command :h3
pair_style multi/lucy/rx style N keyword ... :pre
style = {lookup} or {linear} = method of interpolation
-N = use N values in {lookup}, {linear} tables :ul
+N = use N values in {lookup}, {linear} tables
+weighting = fractional or molecular (optional) :ul
[Examples:]
pair_style multi/lucy/rx linear 1000
+pair_style multi/lucy/rx linear 1000 fractional
+pair_style multi/lucy/rx linear 1000 molecular
pair_coeff * * multibody.table ENTRY1 h2o h2o 7.0
pair_coeff * * multibody.table ENTRY1 h2o 1fluid 7.0 :pre
@@ -94,8 +97,10 @@ tags must either correspond to the species defined in the reaction
kinetics files specified with the "fix rx"_fix_rx.html command or they
must correspond to the tag "1fluid", signifying interaction with a
product species mixture determined through a one-fluid approximation.
-The interaction potential is weighted by the geometric average of the
-concentrations of the two species. The coarse-grained potential is
+The interaction potential is weighted by the geometric average of
+either the mole fraction concentrations or the number of molecules
+associated with the interacting coarse-grained particles (see the
+{fractional} or {molecular} weighting pair style options). The coarse-grained potential is
stored before and after the reaction kinetics solver is applied, where
the difference is defined to be the internal chemical energy (uChem).
@@ -205,7 +210,7 @@ LAMMPS"_Section_start.html#start_3 section for more info.
"pair_coeff"_pair_coeff.html
-[Default:] none
+[Default:] fractional weighting
:line
diff --git a/doc/src/pair_table_rx.txt b/doc/src/pair_table_rx.txt
index e6006f62e22516dd6271a53aa863ecec85cc4e11..d089a4f9da057858b8708a2695f8c45d582177c1 100644
--- a/doc/src/pair_table_rx.txt
+++ b/doc/src/pair_table_rx.txt
@@ -10,16 +10,17 @@ pair_style table/rx command :h3
[Syntax:]
-pair_style table style N :pre
+pair_style table style N ... :pre
style = {lookup} or {linear} or {spline} or {bitmap} = method of interpolation
N = use N values in {lookup}, {linear}, {spline} tables
-N = use 2^N values in {bitmap} tables
+weighting = fractional or molecular (optional) :ul
[Examples:]
pair_style table/rx linear 1000
-pair_style table/rx bitmap 12
+pair_style table/rx linear 1000 fractional
+pair_style table/rx linear 1000 molecular
pair_coeff * * rxn.table ENTRY1 h2o h2o 10.0
pair_coeff * * rxn.table ENTRY1 1fluid 1fluid 10.0
pair_coeff * 3 rxn.table ENTRY1 h2o no2 10.0 :pre
@@ -84,8 +85,10 @@ tags must either correspond to the species defined in the reaction
kinetics files specified with the "fix rx"_fix_rx.html command or they
must correspond to the tag "1fluid", signifying interaction with a
product species mixture determined through a one-fluid approximation.
-The interaction potential is weighted by the geometric average of the
-concentrations of the two species. The coarse-grained potential is
+The interaction potential is weighted by the geometric average of
+either the mole fraction concentrations or the number of molecules
+associated with the interacting coarse-grained particles (see the
+{fractional} or {molecular} weighting pair style options). The coarse-grained potential is
stored before and after the reaction kinetics solver is applied, where
the difference is defined to be the internal chemical energy (uChem).
@@ -230,7 +233,7 @@ LAMMPS"_Section_start.html#start_3 section for more info.
"pair_coeff"_pair_coeff.html
-[Default:] none
+[Default:] fractional weighting
:line
diff --git a/examples/USER/dpd/dpde-vv/log.dpde-vv.reference b/examples/USER/dpd/dpde-vv/log.dpde-vv.reference
index 7bc7bda365b2573c283329239f519f0e14ed143c..800a39f7a548bb7b2e1a2ee20bb0a5c800e34754 100644
--- a/examples/USER/dpd/dpde-vv/log.dpde-vv.reference
+++ b/examples/USER/dpd/dpde-vv/log.dpde-vv.reference
@@ -35,129 +35,133 @@ thermo_modify format float %24.16f
run 1000
Neighbor list info ...
- 1 neighbor list requests
update every 1 steps, delay 0 steps, check no
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 10.6
ghost atom cutoff = 10.6
- binsize = 5.3 -> bins = 25 25 25
-Memory usage per processor = 3.36353 Mbytes
+ binsize = 5.3, bins = 25 25 25
+ 1 neighbor lists, perpetual/occasional/extra = 1 0 0
+ (1) pair dpd/fdt/energy, perpetual
+ pair build: half/bin/newton
+ stencil: half/bin/3d/newton
+ bin: standard
+Memory usage per processor = 4.28221 Mbytes
Step Temp Press PotEng KinEng c_dpdU[1] c_dpdU[2] v_totEnergy c_dpdU[4]
- 0 301.4391322267262012 1636.1776395935085020 1188.6488072196075336 394.4722035796053206 7852.5601874986105031 7852.5601874986105031 17288.2413857964347699 299.9999999999841407
- 10 301.4791572483523510 1486.4422375141198245 1188.7147620806101713 394.5245815119678241 7852.5601874999802021 7852.3731942333779443 17288.1727253259377903 299.9960221120699089
- 20 301.4275643919337426 1677.9356110821624952 1188.7839634625399867 394.4570655673388728 7852.5601874999938445 7852.3711851933012440 17288.1724017231754260 299.9955485734552099
- 30 301.2240988054542186 1452.7304951528931269 1188.8550809767796181 394.1908044563202225 7852.5601875000002110 7852.5679666239848302 17288.1740395570850524 299.9988968405210130
- 40 301.1023506886409677 1527.9758363521380033 1188.9264527568634549 394.0314812537677653 7852.5601874999947540 7852.6574764573806533 17288.1755979680056043 300.0001694462812338
- 50 301.0409654880461972 1597.1737251233498682 1188.9944523606982330 393.9511507566391515 7852.5601875000029395 7852.6700547249911324 17288.1758453423317405 299.9999653064982681
- 60 301.2904978886139133 1610.8630327676828529 1189.0651026961211301 394.2776962691256131 7852.5601874999829306 7852.2734988976435488 17288.1764853628737910 299.9919857290491905
- 70 300.8575037843163500 1489.3259312130880971 1189.1295686642290548 393.7110673208616731 7852.5601874999856591 7852.7707182199101226 17288.1715417049854295 300.0010992278233175
- 80 300.5955830326474825 1449.3896097889587509 1189.1880764967559116 393.3683100440913449 7852.5601875000411383 7853.0484238882281716 17288.1649979291178170 300.0059513551503301
- 90 301.0092332775843147 1553.9266324350364812 1189.2470037925052111 393.9096250433288446 7852.5601875000420478 7852.4452067113825251 17288.1620230472581170 299.9940347326859182
- 100 301.0478004479094238 1539.2270336322194453 1189.3010269201699884 393.9600951881690207 7852.5601875000074870 7852.3416236045995902 17288.1629332129450631 299.9916385566916119
- 110 300.9609384905550087 1500.0429484565006533 1189.3524514939088021 393.8464250502817663 7852.5601874999983920 7852.4114980357189779 17288.1705620799075405 299.9925626482005327
- 120 300.9625536631411933 1630.5065919443034090 1189.4006029528841282 393.8485387131115658 7852.5601875000575092 7852.3600810123671181 17288.1694101784196391 299.9911580775880680
- 130 301.0373750247310340 1539.2267307640183844 1189.4426173625224692 393.9464521696795032 7852.5601874999993015 7852.2178388309775983 17288.1670958631802932 299.9879581026651749
- 140 300.7465104415114752 1550.8353679735087098 1189.4887352231000932 393.5658181350791551 7852.5601874999920256 7852.5559582333216895 17288.1706990914935886 299.9939749909034958
- 150 300.6667173911141617 1634.8987162883277051 1189.5368575067818711 393.4613985788388959 7852.5601874999920256 7852.6079668015609059 17288.1664103871735279 299.9946423938895350
- 160 300.4684731724562425 1462.9400882126803936 1189.5825022927965620 393.2019703048678707 7852.5601874999847496 7852.8265187980177870 17288.1711788956672535 299.9983600613423960
- 170 300.1439323338466920 1510.2352578813552100 1189.6305700279478970 392.7772665220106774 7852.5601874999802021 7853.2009671047335360 17288.1689911546709482 300.0051118582463232
- 180 300.1074244553407198 1529.6307083879951279 1189.6764977580119194 392.7294912276224181 7852.5601874999729262 7853.2047509722533505 17288.1709274578606710 300.0047089238623812
- 190 300.4193298066089142 1546.3205495807171701 1189.7172820166240399 393.1376598363699486 7852.5601874999847496 7852.7461854379371289 17288.1613147909156396 299.9954451643528728
- 200 300.3353919251508728 1532.5496449337254035 1189.7600175880224924 393.0278162310690391 7852.5601874999683787 7852.8107089913455638 17288.1587303104060993 299.9962707550171785
- 210 300.3276568499739483 1504.8178651700843602 1189.7998299597820733 393.0176938818990493 7852.5601875000156724 7852.7810130200659842 17288.1587243617614149 299.9953436245502871
- 220 300.5768315696971626 1592.5896084568344122 1189.8391466344742184 393.3437713226064716 7852.5601875000329528 7852.4205574703573802 17288.1636629274726147 299.9880321846658831
- 230 300.6587445618569063 1672.3049358942289473 1189.8766340798690635 393.4509650976162334 7852.5601874999847496 7852.2733199687863817 17288.1611066462573945 299.9848228571166828
- 240 300.7517707836825025 1527.1722267937811921 1189.9126240081129708 393.5727019751183207 7852.5601875000065775 7852.1160682173085661 17288.1615817005440476 299.9814952182625802
- 250 300.8473715548367409 1589.1847713095248764 1189.9441342461948352 393.6978079843565865 7852.5601875000047585 7851.9625847797888127 17288.1647145103452203 299.9782210858571148
- 260 300.8450266408960942 1623.1896863377055524 1189.9636161513917614 393.6947393603111891 7852.5601874999820211 7851.9471828473988353 17288.1657258590821584 299.9775302202895659
- 270 300.6663619570709898 1564.5160171187899323 1189.9764081239700317 393.4609334472908131 7852.5601875000193104 7852.1708276117251444 17288.1683566830033669 299.9812899253168439
- 280 300.7668534205726019 1618.5400526904263643 1189.9872008155405183 393.5924395618274048 7852.5601875000184009 7852.0271568534708422 17288.1669847308585304 299.9781169783826158
- 290 300.8462727198648849 1562.6765776748122789 1189.9918265985252219 393.6963700162682471 7852.5601875000211294 7851.9189772084127981 17288.1673613232269417 299.9756806168044250
- 300 300.8095414073812890 1525.1785808192844343 1189.9873922767767453 393.6483023295390922 7852.5601875000020300 7851.9657301693578120 17288.1616122756749974 299.9761279889730758
- 310 300.9496330741350221 1566.5597234051326723 1189.9752299662607129 393.8316304464934774 7852.5601875000056680 7851.7898117189633922 17288.1568596317229094 299.9723726900590464
- 320 301.2370566356515837 1513.6869483705047514 1189.9626455872523820 394.2077614578674343 7852.5601874999929350 7851.4248466706330873 17288.1554412157456682 299.9650543775110236
- 330 301.3279721508968692 1549.0667862452519330 1189.9513389477854162 394.3267362020337146 7852.5601874999929350 7851.3129955581916875 17288.1512582080031279 299.9625537201162615
- 340 301.1145736537583844 1414.7930515101759283 1189.9408691169965095 394.0474765890400590 7852.5601874999993015 7851.6028846074832472 17288.1514178135184920 299.9677356565828745
- 350 301.1651600907370039 1529.8016115175887535 1189.9314470205476937 394.1136755032911196 7852.5601874999929350 7851.5441417268757505 17288.1494517507089768 299.9662576716461331
- 360 301.0550563185083206 1536.7721716375504002 1189.9200519814730796 393.9695904359920178 7852.5601875000074870 7851.7101209691463737 17288.1599508866202086 299.9690811750865009
- 370 301.1008976932964742 1522.3385843459479929 1189.9109162496640693 394.0295798208944120 7852.5601875000211294 7851.6603423306560217 17288.1610259012340975 299.9677565060027860
- 380 301.1656898730700505 1505.0548721701993600 1189.9005648244351505 394.1143687921909304 7852.5601875000056680 7851.5816827598300733 17288.1568038764598896 299.9659906785156522
- 390 300.8379322662876802 1740.9151205755624687 1189.8851457594087151 393.6854554509390596 7852.5601875000238579 7852.0268864110385039 17288.1576751214088290 299.9741278188615752
- 400 300.8663790447546376 1564.9461156870302148 1189.8690133470408909 393.7226817503372445 7852.5601875000411383 7852.0043792319993372 17288.1562618294192362 299.9732593416579789
- 410 300.6263441860635908 1564.2840871092373618 1189.8566574093877080 393.4085650033033517 7852.5601874999892971 7852.3284491703725507 17288.1538590830532485 299.9792095875052951
- 420 300.5302259436974168 1438.1569922368764765 1189.8406936554465574 393.2827818158641549 7852.5601875000302243 7852.4696075433648730 17288.1532705147074012 299.9815165752025337
- 430 300.5877786105220935 1503.3641639033023694 1189.8251514530138593 393.3580969454444016 7852.5601874999802021 7852.4023373559457468 17288.1457732543858583 299.9798346272511935
- 440 300.7289160804472772 1689.2527029957295781 1189.8035410609209066 393.5427936314976591 7852.5601875000029395 7852.2436462415198548 17288.1501684339418716 299.9764596782897570
- 450 300.9487198282456575 1497.3668092174791582 1189.7808137689632986 393.8304353457919547 7852.5601874999938445 7851.9788323927432430 17288.1502690074921702 299.9710227473042323
- 460 300.9359942496024587 1625.1573864018491804 1189.7615359247627111 393.8137822755282400 7852.5601875000147629 7852.0165192783370003 17288.1520249786408385 299.9713565393226986
- 470 301.0000133856357252 1486.1561922844011860 1189.7439269526955741 393.8975596188205941 7852.5601874999656502 7851.9561324572268859 17288.1578065287103527 299.9697143418395626
- 480 300.8568627175957886 1535.6080526199095857 1189.7237810071801505 393.7102284019063063 7852.5601874999601932 7852.1697010727630186 17288.1638979818089865 299.9732503057674080
- 490 301.0608040775520067 1497.3221544489886128 1189.7062242497636362 393.9771121242308709 7852.5601874999974825 7851.9258988739011329 17288.1694227478947141 299.9682362511933320
- 500 301.0232592587148019 1517.5854528541199215 1189.6911287485861521 393.9279798589197981 7852.5601875000247674 7851.9823225510326665 17288.1616186585633841 299.9690333355835037
- 510 300.7038579923685120 1420.2615974401142012 1189.6747661513456933 393.5100018730125839 7852.5601874999674692 7852.4114869568047652 17288.1564424811294884 299.9768186576545759
- 520 300.5917863355052759 1537.4862082427132464 1189.6604754398756540 393.3633415734188361 7852.5601875000029395 7852.5789017095057716 17288.1629062228021212 299.9795694302102333
- 530 300.4751352158502868 1481.1071694751799441 1189.6453243069925065 393.2106884527691477 7852.5601874999811116 7852.7451655714066874 17288.1613658311471227 299.9823181268525900
- 540 300.5380123640739498 1547.3461372766389559 1189.6261485232855648 393.2929713568877332 7852.5601875000375003 7852.6850583598352387 17288.1643657400454686 299.9808112190538623
- 550 300.4253885005187499 1544.3485889749692888 1189.6033595464525661 393.1455884232119047 7852.5601874999756546 7852.8598718466746504 17288.1690073163154011 299.9835860164698147
- 560 300.3263552442093101 1556.5150300058251105 1189.5759163336824713 393.0159905619273673 7852.5601875000111249 7853.0148613782675966 17288.1669557738860021 299.9861837797674866
- 570 300.1977324643196425 1511.2320626303917379 1189.5441090918316149 392.8476709710407704 7852.5601875000102154 7853.2098259401755058 17288.1617935030590161 299.9896761688499964
- 580 300.3543631005173893 1588.9566243200433746 1189.5094471319721379 393.0526424747489500 7852.5601875000156724 7853.0374555421631158 17288.1597326488990802 299.9859298211933378
- 590 300.5019108864805730 1504.4406939723214691 1189.4809412920112663 393.2457278908070748 7852.5601874999874781 7852.8704277855340479 17288.1572844683396397 299.9823573257917815
- 600 300.4791158523048011 1540.4690749004150803 1189.4551948503105905 393.2158976318902432 7852.5601875000220389 7852.9312239063838206 17288.1625038886049879 299.9832002920041987
- 610 300.5939139841889869 1368.0565839211087678 1189.4252547652590692 393.3661258776944578 7852.5601874999574648 7852.8130977336286378 17288.1646658765384927 299.9807742697515778
- 620 300.7674247480806002 1483.2566452708945235 1189.3941250938435132 393.5931872179773450 7852.5601875000193104 7852.6187967208716145 17288.1662965327122947 299.9766963671718258
- 630 300.7920034341021278 1543.0699124130637756 1189.3598279316649950 393.6253516166882491 7852.5601875000302243 7852.6219971866230480 17288.1673642350069713 299.9762538437230432
- 640 300.8032734267029014 1423.2549819291616586 1189.3293074476885067 393.6400998638143278 7852.5601874999847496 7852.6384826097782934 17288.1680774212654796 299.9762118202994543
- 650 300.7516995878241346 1542.6559695158523482 1189.3021161045705867 393.5726088061030055 7852.5601874999720167 7852.7361949473242930 17288.1711073579681397 299.9775656396505497
- 660 300.8699697098109596 1675.5121937767839881 1189.2687179804190691 393.7273806013013768 7852.5601874999802021 7852.6179739687149777 17288.1742600504148868 299.9750492262036801
- 670 301.0255004186900578 1520.7397686587873977 1189.2284265783687260 393.9309127074437242 7852.5601874999847496 7852.4592279727157802 17288.1787547585117863 299.9715123049731460
- 680 301.1071983488760679 1651.9751417063259851 1189.1858967311386550 394.0378250459656329 7852.5601875000002110 7852.3982826328638112 17288.1821919099675142 299.9699481289110850
- 690 301.0027086454253435 1496.1607274163641250 1189.1436949551202815 393.9010867158519886 7852.5601875000293148 7852.5788938360938118 17288.1838630070960789 299.9731939774295597
- 700 300.9009090279179759 1551.8182127127668082 1189.0993919251338866 393.7678687121208441 7852.5601875000102154 7852.7513665452252098 17288.1788146824910655 299.9761043445071209
- 710 301.2325536720837817 1678.1546953970853338 1189.0528341066981284 394.2018687459686817 7852.5601874999956635 7852.3633298995819132 17288.1782202522445004 299.9683013583347133
- 720 301.2122298224125529 1524.1415452491430642 1189.0046957644285612 394.1752723525083866 7852.5601875000093059 7852.4351629896145823 17288.1753186065616319 299.9693315350040734
- 730 301.0763282392692304 1547.1987029633166912 1188.9602551214045434 393.9974275034455218 7852.5601874999883876 7852.6518053705112834 17288.1696754953518393 299.9732715774841267
- 740 301.3262401480515109 1544.7045314021493141 1188.9131307177485724 394.3244696516559884 7852.5601874999965730 7852.3694201272974169 17288.1672079966992897 299.9674666811455950
- 750 301.5740779122830304 1591.1785078054851965 1188.8637580645938669 394.6487975126887022 7852.5601875000029395 7852.0919529470393172 17288.1646960243233480 299.9616008527094095
- 760 301.4385361878654521 1547.3218422039201414 1188.8113669183098864 394.4714235854450521 7852.5601874999838401 7852.3161911124070684 17288.1591691161447670 299.9656339783694534
- 770 301.6110125684814420 1494.5039561806622714 1188.7581685915934031 394.6971313010439530 7852.5601875000083965 7852.1351720579104949 17288.1506594505553949 299.9619855799395509
- 780 301.8360352039435384 1588.1458619705292676 1188.7039178696472845 394.9916026067776329 7852.5601874999956635 7851.9015195838428554 17288.1572275602629816 299.9572350302977952
- 790 302.1008324754310479 1545.4409171812178556 1188.6491103416560691 395.3381241828382144 7852.5601875000138534 7851.6150048936624444 17288.1624269181702402 299.9513959104631340
- 800 301.9660372380565718 1563.9565804790736365 1188.5964649891604950 395.1617271307158035 7852.5601874999874781 7851.8461249560614306 17288.1645045759250934 299.9555810527747326
- 810 302.0507207347627627 1511.4560763489957935 1188.5468477146612258 395.2725464702810996 7852.5601875000120344 7851.7904104899025697 17288.1699921748586348 299.9541551776504775
- 820 302.4700213214911741 1458.5135514273570152 1188.4981381693974072 395.8212556746473751 7852.5601875000202199 7851.2935886962204677 17288.1731700402851857 299.9441803241180651
- 830 302.2853997979337350 1496.2544527963129894 1188.4496917372191547 395.5796544641875698 7852.5601875000447762 7851.5862641793482908 17288.1757978808018379 299.9494768794835977
- 840 302.0840465730901201 1518.8301331998704882 1188.3994383226176978 395.3161576523596636 7852.5601875000038490 7851.8962146812327774 17288.1719981562127941 299.9550476592922337
- 850 301.8910942560261788 1469.8827850510901953 1188.3489956121345585 395.0636545180261692 7852.5601874999829306 7852.2025804631493884 17288.1754180932912277 299.9606927700139067
- 860 301.7284384160519153 1657.6802015862324424 1188.3052233777652873 394.8507982536594341 7852.5601875000093059 7852.4644669022691232 17288.1806760337058222 299.9652835238809985
- 870 301.6331619894115192 1501.5829953208524330 1188.2628815714097072 394.7261166912876433 7852.5601875000202199 7852.6378180648598573 17288.1870038275774277 299.9682811831179379
- 880 301.3703918424367316 1499.1595903074553462 1188.2195190931643083 394.3822478705861272 7852.5601874999956635 7853.0266423250832304 17288.1885967888301820 299.9755099056966401
- 890 301.4157954313303662 1598.8758859042511631 1188.1845892608291706 394.4416643558612918 7852.5601875000065775 7853.0036606192506952 17288.1901017359487014 299.9745322513492738
- 900 301.4752150615485675 1621.2148728756822038 1188.1517520946135846 394.5194226492019993 7852.5601874999711072 7852.9579580608560718 17288.1893203046420240 299.9733125337182287
- 910 301.4308816315938770 1538.4823217911632582 1188.1159856659232901 394.4614066057066566 7852.5601875000002110 7853.0558695713261841 17288.1934493429580471 299.9748317405193916
- 920 301.4323110133492492 1594.7193046491217956 1188.0835779842032025 394.4632771371357762 7852.5601875000202199 7853.0942701464364291 17288.2013127677964803 299.9751127806911200
- 930 301.4801256941950101 1387.6885377097617038 1188.0464206196895702 394.5258488489681099 7852.5601875000229484 7853.0656502842994087 17288.1981072529815719 299.9740698440909910
- 940 301.8075611840245074 1534.2487040663793323 1188.0124217312886685 394.9543406584059539 7852.5601874999701977 7852.6729444202819650 17288.1998943099461030 299.9660570413493588
- 950 301.6915970126173647 1567.7725992489238251 1187.9790455470049437 394.8025864986412898 7852.5601875000274958 7852.8619557087595240 17288.2037752544347313 299.9694678653150959
- 960 301.6392594677008105 1504.8502165144939227 1187.9439133338105421 394.7340960325207675 7852.5601874999711072 7852.9728807988849439 17288.2110776651898050 299.9711546356286362
- 970 301.6049535791644303 1514.0198965433548892 1187.9094123369413865 394.6892023276233772 7852.5601874999765641 7853.0497909819878259 17288.2085931465298927 299.9722547114341751
- 980 301.2982841679705643 1634.1208149125807267 1187.8768454876480973 394.2878856256063500 7852.5601874999856591 7853.4862008383515786 17288.2111194515891839 299.9802110109069986
- 990 301.2573007350166563 1489.7316698898257528 1187.8432331161868660 394.2342534877078606 7852.5601875000047585 7853.5840096862748396 17288.2216837901723920 299.9819468620868292
- 1000 301.3195135766228532 1562.6587211933920116 1187.8034267774903583 394.3156670604516307 7852.5601874999356369 7853.5372636956635688 17288.2165450335414789 299.9807651637231629
-Loop time of 21.3308 on 1 procs for 1000 steps with 10125 atoms
+ 0 301.4391322267262012 1636.1776395935080473 1188.6488072196075336 394.4722035796053206 0.0000000000000000 15705.1203749972210062 17288.2413857964347699 299.9999999999841407
+ 10 301.4791572483523510 1486.4422375141214161 1188.7147620806101713 394.5245815119678241 0.0000000000000000 15704.9333817333845218 17288.1727253259632562 299.9960221120699089
+ 20 301.4275643919337995 1677.9356110821622678 1188.7839634625399867 394.4570655673389865 -0.0000000000000000 15704.9313726932996360 17288.1724017231790640 299.9955485734552667
+ 30 301.2240988054542186 1452.7304951528922174 1188.8550809767796181 394.1908044563202225 -0.0000000000000000 15705.1281541239713988 17288.1740395570705005 299.9988968405209562
+ 40 301.1023506886409109 1527.9758363521384581 1188.9264527568634549 394.0314812537677085 -0.0000000000000000 15705.2176639573335706 17288.1755979679655866 300.0001694462812907
+ 50 301.0409654880461972 1597.1737251233505503 1188.9944523606984603 393.9511507566391515 -0.0000000000000000 15705.2302422249904339 17288.1758453423281026 299.9999653064982112
+ 60 301.2904978886138565 1610.8630327676828529 1189.0651026961211301 394.2776962691255562 -0.0000000000000000 15704.8336863976528548 17288.1764853628992569 299.9919857290491905
+ 70 300.8575037843164068 1489.3259312130892340 1189.1295686642290548 393.7110673208617300 0.0000000000000000 15705.3309057198275696 17288.1715417049199459 300.0010992278232607
+ 80 300.5955830326474825 1449.3896097889576140 1189.1880764967559116 393.3683100440913449 -0.0000000000000000 15705.6086113882302016 17288.1649979290777992 300.0059513551502164
+ 90 301.0092332775843147 1553.9266324350371633 1189.2470037925056658 393.9096250433288446 -0.0000000000000000 15705.0053942113881931 17288.1620230472217372 299.9940347326859182
+ 100 301.0478004479094238 1539.2270336322201274 1189.3010269201699884 393.9600951881690207 -0.0000000000000000 15704.9018111045588739 17288.1629332128977694 299.9916385566916119
+ 110 300.9609384905550655 1500.0429484565015628 1189.3524514939088021 393.8464250502818231 -0.0000000000000000 15704.9716855356964516 17288.1705620798857126 299.9925626482006464
+ 120 300.9625536631413070 1630.5065919443020448 1189.4006029528841282 393.8485387131116795 0.0000000000000000 15704.9202685123345873 17288.1694101783286897 299.9911580775880680
+ 130 301.0373750247309772 1539.2267307640188392 1189.4426173625224692 393.9464521696794463 -0.0000000000000000 15704.7780263310032751 17288.1670958632057591 299.9879581026650044
+ 140 300.7465104415114183 1550.8353679735089372 1189.4887352231000932 393.5658181350790983 0.0000000000000000 15705.1161457332873397 17288.1706990914681228 299.9939749909034958
+ 150 300.6667173911142186 1634.8987162883267956 1189.5368575067818711 393.4613985788390096 0.0000000000000000 15705.1681543015274656 17288.1664103871480620 299.9946423938894213
+ 160 300.4684731724561857 1462.9400882126797114 1189.5825022927965620 393.2019703048678139 0.0000000000000000 15705.3867062980680203 17288.1711788957327371 299.9983600613422254
+ 170 300.1439323338466920 1510.2352578813547552 1189.6305700279476696 392.7772665220106774 -0.0000000000000000 15705.7611546046609874 17288.1689911546200165 300.0051118582463232
+ 180 300.1074244553407766 1529.6307083879964921 1189.6764977580119194 392.7294912276225318 -0.0000000000000000 15705.7649384723172261 17288.1709274579516205 300.0047089238623812
+ 190 300.4193298066088573 1546.3205495807169427 1189.7172820166242673 393.1376598363698349 0.0000000000000000 15705.3063729379555298 17288.1613147909483814 299.9954451643527022
+ 200 300.3353919251508728 1532.5496449337249487 1189.7600175880224924 393.0278162310690391 -0.0000000000000000 15705.3708964914076205 17288.1587303105006868 299.9962707550172922
+ 210 300.3276568499739483 1504.8178651700850423 1189.7998299597820733 393.0176938818990493 0.0000000000000000 15705.3412005200552812 17288.1587243617359491 299.9953436245502871
+ 220 300.5768315696972195 1592.5896084568353217 1189.8391466344739911 393.3437713226065284 -0.0000000000000000 15704.9807449702821032 17288.1636629273634753 299.9880321846658262
+ 230 300.6587445618569063 1672.3049358942282652 1189.8766340798690635 393.4509650976162334 0.0000000000000000 15704.8335074687693123 17288.1611066462537565 299.9848228571169102
+ 240 300.7517707836825025 1527.1722267937814195 1189.9126240081131982 393.5727019751183207 -0.0000000000000000 15704.6762557172896777 17288.1615817005222198 299.9814952182625802
+ 250 300.8473715548367409 1589.1847713095232848 1189.9441342461948352 393.6978079843565865 0.0000000000000000 15704.5227722798481409 17288.1647145103997900 299.9782210858571148
+ 260 300.8450266408959806 1623.1896863377055524 1189.9636161513917614 393.6947393603110186 0.0000000000000000 15704.5073703474117792 17288.1657258591149002 299.9775302202894522
+ 270 300.6663619570710466 1564.5160171187892502 1189.9764081239700317 393.4609334472908699 0.0000000000000000 15704.7310151116998895 17288.1683566829597112 299.9812899253167302
+ 280 300.7668534205727155 1618.5400526904256822 1189.9872008155405183 393.5924395618275184 0.0000000000000000 15704.5873443533891987 17288.1669847307566670 299.9781169783825590
+ 290 300.8462727198648281 1562.6765776748138705 1189.9918265985252219 393.6963700162681334 0.0000000000000000 15704.4791647084566648 17288.1673613232487696 299.9756806168042544
+ 300 300.8095414073812890 1525.1785808192844343 1189.9873922767767453 393.6483023295390922 0.0000000000000000 15704.5259176693853078 17288.1616122757004632 299.9761279889731327
+ 310 300.9496330741349652 1566.5597234051326723 1189.9752299662607129 393.8316304464933637 0.0000000000000000 15704.3499992189717887 17288.1568596317265474 299.9723726900589327
+ 320 301.2370566356514132 1513.6869483705036146 1189.9626455872523820 394.2077614578672069 0.0000000000000000 15703.9850341706151085 17288.1554412157347542 299.9650543775107394
+ 330 301.3279721508969260 1549.0667862452526151 1189.9513389477854162 394.3267362020338282 0.0000000000000000 15703.8731830581982649 17288.1512582080176799 299.9625537201162615
+ 340 301.1145736537582707 1414.7930515101757010 1189.9408691169962822 394.0474765890398885 0.0000000000000000 15704.1630721074998291 17288.1514178135366819 299.9677356565827040
+ 350 301.1651600907369470 1529.8016115175894356 1189.9314470205474663 394.1136755032910628 0.0000000000000000 15704.1043292268568621 17288.1494517506944248 299.9662576716459625
+ 360 301.0550563185083206 1536.7721716375513097 1189.9200519814730796 393.9695904359920178 0.0000000000000000 15704.2703084691693221 17288.1599508866347605 299.9690811750866146
+ 370 301.1008976932965311 1522.3385843459491298 1189.9109162496640693 394.0295798208944689 0.0000000000000000 15704.2205298306434997 17288.1610259012013557 299.9677565060027860
+ 380 301.1656898730701073 1505.0548721701995873 1189.9005648244356053 394.1143687921909873 -0.0000000000000000 15704.1418702597857191 17288.1568038764125959 299.9659906785157091
+ 390 300.8379322662877371 1740.9151205755633782 1189.8851457594089425 393.6854554509391164 -0.0000000000000000 15704.5870739109432179 17288.1576751212924137 299.9741278188614046
+ 400 300.8663790447545239 1564.9461156870302148 1189.8690133470406636 393.7226817503371308 0.0000000000000000 15704.5645667319495260 17288.1562618293282867 299.9732593416576947
+ 410 300.6263441860637045 1564.2840871092375892 1189.8566574093874806 393.4085650033035222 -0.0000000000000000 15704.8886366703736712 17288.1538590830641624 299.9792095875053519
+ 420 300.5302259436973031 1438.1569922368769312 1189.8406936554461026 393.2827818158640412 0.0000000000000000 15705.0297950433650840 17288.1532705146746594 299.9815165752024768
+ 430 300.5877786105221503 1503.3641639033021420 1189.8251514530136319 393.3580969454445153 -0.0000000000000000 15704.9625248558968451 17288.1457732543567545 299.9798346272512504
+ 440 300.7289160804472772 1689.2527029957295781 1189.8035410609209066 393.5427936314976591 -0.0000000000000000 15704.8038337415237038 17288.1501684339418716 299.9764596782894728
+ 450 300.9487198282456006 1497.3668092174784761 1189.7808137689632986 393.8304353457918978 -0.0000000000000000 15704.5390198927143501 17288.1502690074703423 299.9710227473042323
+ 460 300.9359942496024019 1625.1573864018473614 1189.7615359247631659 393.8137822755281263 0.0000000000000000 15704.5767067783035600 17288.1520249785935448 299.9713565393225849
+ 470 301.0000133856357252 1486.1561922844020955 1189.7439269526958014 393.8975596188205941 0.0000000000000000 15704.5163199572089070 17288.1578065287249046 299.9697143418395058
+ 480 300.8568627175958454 1535.6080526199100404 1189.7237810071803779 393.7102284019064200 -0.0000000000000000 15704.7298885727686866 17288.1638979818562802 299.9732503057675785
+ 490 301.0608040775520067 1497.3221544489890675 1189.7062242497640909 393.9771121242308709 -0.0000000000000000 15704.4860863739140768 17288.1694227479092660 299.9682362511933889
+ 500 301.0232592587148019 1517.5854528541185573 1189.6911287485863795 393.9279798589197981 -0.0000000000000000 15704.5425100510510674 17288.1616186585561081 299.9690333355832195
+ 510 300.7038579923685120 1420.2615974401142012 1189.6747661513456933 393.5100018730125839 -0.0000000000000000 15704.9716744568013382 17288.1564424811585923 299.9768186576548032
+ 520 300.5917863355052759 1537.4862082427125642 1189.6604754398761088 393.3633415734188361 -0.0000000000000000 15705.1390892093895673 17288.1629062226857059 299.9795694302102902
+ 530 300.4751352158504574 1481.1071694751785799 1189.6453243069920518 393.2106884527693751 -0.0000000000000000 15705.3053530714041699 17288.1613658311653126 299.9823181268525900
+ 540 300.5380123640739498 1547.3461372766387285 1189.6261485232855648 393.2929713568877332 0.0000000000000000 15705.2452458598490921 17288.1643657400236407 299.9808112190538623
+ 550 300.4253885005187499 1544.3485889749688340 1189.6033595464525661 393.1455884232119047 0.0000000000000000 15705.4200593467012368 17288.1690073163663328 299.9835860164698147
+ 560 300.3263552442091395 1556.5150300058239736 1189.5759163336820166 393.0159905619271399 0.0000000000000000 15705.5750488783432957 17288.1669557739514858 299.9861837797674298
+ 570 300.1977324643196994 1511.2320626303924200 1189.5441090918316149 392.8476709710408272 0.0000000000000000 15705.7700134401693504 17288.1617935030408262 299.9896761688500533
+ 580 300.3543631005173893 1588.9566243200420104 1189.5094471319723652 393.0526424747489500 -0.0000000000000000 15705.5976430422142585 17288.1597326489354600 299.9859298211932810
+ 590 300.5019108864805730 1504.4406939723210144 1189.4809412920112663 393.2457278908070748 -0.0000000000000000 15705.4306152855297114 17288.1572844683469157 299.9823573257918952
+ 600 300.4791158523048011 1540.4690749004137160 1189.4551948503108179 393.2158976318902432 0.0000000000000000 15705.4914114063831221 17288.1625038885831600 299.9832002920041418
+ 610 300.5939139841890437 1368.0565839211083130 1189.4252547652597514 393.3661258776945715 0.0000000000000000 15705.3732852337052464 17288.1646658766585460 299.9807742697515209
+ 620 300.7674247480806002 1483.2566452708929319 1189.3941250938437406 393.5931872179773450 0.0000000000000000 15705.1789842209145718 17288.1662965327341226 299.9766963671719395
+ 630 300.7920034341022415 1543.0699124130630935 1189.3598279316649950 393.6253516166883628 -0.0000000000000000 15705.1821846865786938 17288.1673642349305737 299.9762538437231001
+ 640 300.8032734267029014 1423.2549819291609765 1189.3293074476887341 393.6400998638143278 -0.0000000000000000 15705.1986701098048798 17288.1680774213091354 299.9762118202993975
+ 650 300.7516995878240209 1542.6559695158514387 1189.3021161045703593 393.5726088061028349 0.0000000000000000 15705.2963824473390559 17288.1711073580117954 299.9775656396504360
+ 660 300.8699697098108459 1675.5121937767842155 1189.2687179804192965 393.7273806013012063 0.0000000000000000 15705.1781614686860848 17288.1742600504076108 299.9750492262035095
+ 670 301.0255004186899441 1520.7397686587889893 1189.2284265783694082 393.9309127074436105 0.0000000000000000 15705.0194154727287241 17288.1787547585408902 299.9715123049731460
+ 680 301.1071983488761248 1651.9751417063253029 1189.1858967311388824 394.0378250459656897 0.0000000000000000 15704.9584701329349627 17288.1821919100402738 299.9699481289110281
+ 690 301.0027086454255141 1496.1607274163641250 1189.1436949551202815 393.9010867158522160 0.0000000000000000 15705.1390813360922039 17288.1838630070633371 299.9731939774292755
+ 700 300.9009090279178622 1551.8182127127668082 1189.0993919251338866 393.7678687121206735 -0.0000000000000000 15705.3115540452217829 17288.1788146824765136 299.9761043445070641
+ 710 301.2325536720837817 1678.1546953970841969 1189.0528341066981284 394.2018687459686817 0.0000000000000000 15704.9235173995584773 17288.1782202522263105 299.9683013583346565
+ 720 301.2122298224125529 1524.1415452491437463 1189.0046957644283339 394.1752723525083866 0.0000000000000000 15704.9953504895402148 17288.1753186064779584 299.9693315350040734
+ 730 301.0763282392692304 1547.1987029633176007 1188.9602551214045434 393.9974275034455218 0.0000000000000000 15705.2119928705469647 17288.1696754953954951 299.9732715774840699
+ 740 301.3262401480515109 1544.7045314021493141 1188.9131307177485724 394.3244696516559884 0.0000000000000000 15704.9296076272603386 17288.1672079966665478 299.9674666811455950
+ 750 301.5740779122830872 1591.1785078054849691 1188.8637580645940943 394.6487975126887591 0.0000000000000000 15704.6521404470349808 17288.1646960243160720 299.9616008527092959
+ 760 301.4385361878655658 1547.3218422039212783 1188.8113669183098864 394.4714235854451658 0.0000000000000000 15704.8763786124927719 17288.1591691162466304 299.9656339783693966
+ 770 301.6110125684815557 1494.5039561806624988 1188.7581685915934031 394.6971313010441236 0.0000000000000000 15704.6953595579507237 17288.1506594505881367 299.9619855799396646
+ 780 301.8360352039435384 1588.1458619705304045 1188.7039178696477393 394.9916026067776329 0.0000000000000000 15704.4617070838321524 17288.1572275602593436 299.9572350302976247
+ 790 302.1008324754310479 1545.4409171812180830 1188.6491103416560691 395.3381241828382144 0.0000000000000000 15704.1751923936917592 17288.1624269181847922 299.9513959104630771
+ 800 301.9660372380565718 1563.9565804790738639 1188.5964649891604950 395.1617271307158035 0.0000000000000000 15704.4063124560707365 17288.1645045759469212 299.9555810527747326
+ 810 302.0507207347627059 1511.4560763489960209 1188.5468477146607711 395.2725464702810427 0.0000000000000000 15704.3505979898400255 17288.1699921747822373 299.9541551776507617
+ 820 302.4700213214913447 1458.5135514273563331 1188.4981381693974072 395.8212556746476025 0.0000000000000000 15703.8537761962070363 17288.1731700402524439 299.9441803241177809
+ 830 302.2853997979336214 1496.2544527963145811 1188.4496917372191547 395.5796544641873993 0.0000000000000000 15704.1464516793694202 17288.1757978807763720 299.9494768794834840
+ 840 302.0840465730901201 1518.8301331998702608 1188.3994383226179252 395.3161576523596636 0.0000000000000000 15704.4564021812439023 17288.1719981562200701 299.9550476592922337
+ 850 301.8910942560260082 1469.8827850510904227 1188.3489956121347859 395.0636545180259986 0.0000000000000000 15704.7627679631386854 17288.1754180932985037 299.9606927700136794
+ 860 301.7284384160518016 1657.6802015862315329 1188.3052233777652873 394.8507982536592635 0.0000000000000000 15705.0246544022065791 17288.1806760336330626 299.9652835238807711
+ 870 301.6331619894114624 1501.5829953208508414 1188.2628815714099346 394.7261166912875865 0.0000000000000000 15705.1980055648327834 17288.1870038275301340 299.9682811831179947
+ 880 301.3703918424367316 1499.1595903074555736 1188.2195190931643083 394.3822478705861272 0.0000000000000000 15705.5868298250898079 17288.1885967888410960 299.9755099056964127
+ 890 301.4157954313303662 1598.8758859042509357 1188.1845892608291706 394.4416643558612918 0.0000000000000000 15705.5638481192290783 17288.1901017359195976 299.9745322513492738
+ 900 301.4752150615486812 1621.2148728756842502 1188.1517520946144941 394.5194226492021699 0.0000000000000000 15705.5181455608308170 17288.1893203046492999 299.9733125337182287
+ 910 301.4308816315937634 1538.4823217911621214 1188.1159856659228353 394.4614066057064861 0.0000000000000000 15705.6160570713091147 17288.1934493429398572 299.9748317405192779
+ 920 301.4323110133492492 1594.7193046491240693 1188.0835779842032025 394.4632771371357762 0.0000000000000000 15705.6544576464475540 17288.2013127677855664 299.9751127806913473
+ 930 301.4801256941949532 1387.6885377097596574 1188.0464206196900250 394.5258488489680531 0.0000000000000000 15705.6258377843460039 17288.1981072530033998 299.9740698440912183
+ 940 301.8075611840245074 1534.2487040663797870 1188.0124217312888959 394.9543406584059539 0.0000000000000000 15705.2331319202457962 17288.1998943099388271 299.9660570413491882
+ 950 301.6915970126175353 1567.7725992489226883 1187.9790455470049437 394.8025864986415172 0.0000000000000000 15705.4221432087451831 17288.2037752543910756 299.9694678653152096
+ 960 301.6392594677008105 1504.8502165144939227 1187.9439133338107695 394.7340960325207675 0.0000000000000000 15705.5330682989206252 17288.2110776652516506 299.9711546356285226
+ 970 301.6049535791644871 1514.0198965433535250 1187.9094123369409317 394.6892023276234909 0.0000000000000000 15705.6099784820144123 17288.2085931465771864 299.9722547114341751
+ 980 301.2982841679706780 1634.1208149125800446 1187.8768454876478700 394.2878856256065205 0.0000000000000000 15706.0463883383199573 17288.2111194515746320 299.9802110109068849
+ 990 301.2573007350166563 1489.7316698898262075 1187.8432331161866387 394.2342534877078606 0.0000000000000000 15706.1441971863041545 17288.2216837901978579 299.9819468620868292
+ 1000 301.3195135766228532 1562.6587211933931485 1187.8034267774903583 394.3156670604516307 0.0000000000000000 15706.0974511956701463 17288.2165450336106005 299.9807651637235040
+Loop time of 17.0881 on 1 procs for 1000 steps with 10125 atoms
-Performance: 4.050 ns/day, 5.925 hours/ns, 46.880 timesteps/s
-99.8% CPU use with 1 MPI tasks x no OpenMP threads
+Performance: 5.056 ns/day, 4.747 hours/ns, 58.520 timesteps/s
+100.0% 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.099 | 10.099 | 10.099 | 0.0 | 47.34
-Neigh | 10.145 | 10.145 | 10.145 | 0.0 | 47.56
-Comm | 0.49807 | 0.49807 | 0.49807 | 0.0 | 2.33
-Output | 0.011203 | 0.011203 | 0.011203 | 0.0 | 0.05
-Modify | 0.28296 | 0.28296 | 0.28296 | 0.0 | 1.33
-Other | | 0.295 | | | 1.38
+Pair | 8.0541 | 8.0541 | 8.0541 | 0.0 | 47.13
+Neigh | 8.1306 | 8.1306 | 8.1306 | 0.0 | 47.58
+Comm | 0.39415 | 0.39415 | 0.39415 | 0.0 | 2.31
+Output | 0.01103 | 0.01103 | 0.01103 | 0.0 | 0.06
+Modify | 0.24061 | 0.24061 | 0.24061 | 0.0 | 1.41
+Other | | 0.2576 | | | 1.51
Nlocal: 10125 ave 10125 max 10125 min
Histogram: 1 0 0 0 0 0 0 0 0 0
@@ -170,4 +174,4 @@ Total # of neighbors = 114682
Ave neighs/atom = 11.3266
Neighbor list builds = 1000
Dangerous builds not checked
-Total wall time: 0:00:21
+Total wall time: 0:00:17
diff --git a/examples/USER/dpd/dpdrx-shardlow/in.dpdrx-shardlow b/examples/USER/dpd/dpdrx-shardlow/in.dpdrx-shardlow
index e65b5a14dbf32a4a30aaebedf81d13474be20c74..815c974741d674fc74af5a025fe5b598ba93894e 100755
--- a/examples/USER/dpd/dpdrx-shardlow/in.dpdrx-shardlow
+++ b/examples/USER/dpd/dpdrx-shardlow/in.dpdrx-shardlow
@@ -37,7 +37,7 @@ timestep 0.001
pair_style hybrid/overlay dpd/fdt/energy 16.00 234324 exp6/rx 16.00
pair_coeff * * dpd/fdt/energy 0.0 0.05 10.0 16.00
-pair_coeff * * exp6/rx params.exp6 1fluid 1fluid 1.0 1.0 16.00
+pair_coeff * * exp6/rx params.exp6 1fluid 1fluid exponent 1.0 1.0 16.00
fix 1 all shardlow
fix 2 all nve
diff --git a/examples/USER/dpd/dpdrx-shardlow/log.dpdrx-shardlow.reference b/examples/USER/dpd/dpdrx-shardlow/log.dpdrx-shardlow.reference
index 067708154a012c642fcd05d0ee94813ef43375cd..b80e033eb9d15d94faca04c0d1ad490b3a3cd372 100644
--- a/examples/USER/dpd/dpdrx-shardlow/log.dpdrx-shardlow.reference
+++ b/examples/USER/dpd/dpdrx-shardlow/log.dpdrx-shardlow.reference
@@ -48,7 +48,7 @@ timestep 0.001
pair_style hybrid/overlay dpd/fdt/energy 16.00 234324 exp6/rx 16.00
pair_coeff * * dpd/fdt/energy 0.0 0.05 10.0 16.00
-pair_coeff * * exp6/rx params.exp6 1fluid 1fluid 1.0 1.0 16.00
+pair_coeff * * exp6/rx params.exp6 1fluid 1fluid exponent 1.0 1.0 16.00
fix 1 all shardlow
fix 2 all nve
@@ -69,39 +69,51 @@ dump_modify 2 sort id
run 10
Neighbor list info ...
- 2 neighbor list requests
update every 1 steps, delay 10 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 18
ghost atom cutoff = 18
- binsize = 9 -> bins = 8 8 5
-Memory usage per processor = 6.52436 Mbytes
+ binsize = 9, bins = 8 8 5
+ 3 neighbor lists, perpetual/occasional/extra = 3 0 0
+ (1) pair dpd/fdt/energy, perpetual
+ pair build: half/bin/newton
+ stencil: half/bin/3d/newton
+ bin: standard
+ (2) pair exp6/rx, perpetual, copy from (1)
+ pair build: copy
+ stencil: none
+ bin: none
+ (3) fix shardlow, perpetual, ssa
+ pair build: half/bin/newton/ssa
+ stencil: half/bin/3d/newton/ssa
+ bin: ssa
+Memory usage per processor = 8.39564 Mbytes
Step Temp Press Volume PotEng KinEng c_dpdU[1] c_dpdU[2] c_dpdU[3] v_totEnergy c_dpdU[4]
- 0 2065.00000000 1368.17463335 179834.51777865 0.00000000 230.35385810 3841.42393279 3841.42393279 0.00000000 7682.84786557 2065.00000000
- 1 2064.93210437 1368.12964881 179834.51777865 0.00000000 230.34628424 3841.42393279 3841.43150665 0.00000000 7682.85543943 2065.20275230
- 2 2067.82089565 1370.04362990 179834.51777865 -0.00000000 230.66853326 3841.42393279 3841.10925763 0.00000000 7682.53319042 2065.32453473
- 3 2070.45225169 1371.78704616 179834.51777865 -0.00000000 230.96206499 3841.42393279 3840.81572590 0.00000000 7682.23965869 2065.45336917
- 4 2075.00241157 1374.80177416 179834.51777865 -0.00000000 231.46964217 3841.42393279 3840.30814872 0.00000000 7681.73208151 2065.52973333
- 5 2073.96509212 1374.11449370 179834.51777865 -0.00000000 231.35392762 3841.42393279 3840.42386327 0.00000000 7681.84779605 2065.76011517
- 6 2074.26516936 1374.31331117 179834.51777865 -0.00000000 231.38740169 3841.42393279 3840.39038920 0.00000000 7681.81432198 2065.95399323
- 7 2071.41069700 1372.42206822 179834.51777865 -0.00000000 231.06898100 3841.42393279 3840.70880989 0.00000000 7682.13274267 2066.23407076
- 8 2071.35844957 1372.38745146 179834.51777865 -0.00000000 231.06315272 3841.42393279 3840.71463817 0.00000000 7682.13857095 2066.43766287
- 9 2071.35676496 1372.38633532 179834.51777865 -0.00000000 231.06296480 3841.42393279 3840.71482609 0.00000000 7682.13875887 2066.64001166
- 10 2066.53172340 1369.18948328 179834.51777865 -0.00000000 230.52472415 3841.42393279 3841.25306673 0.00000000 7682.67699952 2066.97516855
-Loop time of 0.289778 on 1 procs for 10 steps with 864 atoms
-
-Performance: 2.982 ns/day, 8.049 hours/ns, 34.509 timesteps/s
-99.4% CPU use with 1 MPI tasks x no OpenMP threads
+ 0 2065.00000000 1368.17463335 179834.51777865 0.00000000 230.35385810 0.00000000 7682.84786557 0.00000000 7682.84786557 2065.00000000
+ 1 2064.93210437 1368.12964881 179834.51777865 0.00000000 230.34628424 0.00000000 7682.85543943 0.00000000 7682.85543943 2065.20275230
+ 2 2067.82089565 1370.04362990 179834.51777865 -0.00000000 230.66853326 0.00000000 7682.53319042 0.00000000 7682.53319042 2065.32453473
+ 3 2070.45225169 1371.78704616 179834.51777865 -0.00000000 230.96206499 0.00000000 7682.23965869 0.00000000 7682.23965869 2065.45336917
+ 4 2075.00241157 1374.80177416 179834.51777865 -0.00000000 231.46964217 0.00000000 7681.73208151 0.00000000 7681.73208151 2065.52973333
+ 5 2073.96509212 1374.11449370 179834.51777865 -0.00000000 231.35392762 -0.00000000 7681.84779605 0.00000000 7681.84779605 2065.76011517
+ 6 2074.26516936 1374.31331117 179834.51777865 -0.00000000 231.38740169 -0.00000000 7681.81432198 0.00000000 7681.81432198 2065.95399323
+ 7 2071.41069700 1372.42206822 179834.51777865 -0.00000000 231.06898100 -0.00000000 7682.13274267 0.00000000 7682.13274267 2066.23407076
+ 8 2071.35844957 1372.38745146 179834.51777865 -0.00000000 231.06315272 0.00000000 7682.13857095 0.00000000 7682.13857095 2066.43766287
+ 9 2071.35676496 1372.38633532 179834.51777865 -0.00000000 231.06296480 0.00000000 7682.13875887 0.00000000 7682.13875887 2066.64001166
+ 10 2066.53172340 1369.18948328 179834.51777865 -0.00000000 230.52472415 0.00000000 7682.67699952 0.00000000 7682.67699952 2066.97516855
+Loop time of 0.611304 on 1 procs for 10 steps with 864 atoms
+
+Performance: 1.413 ns/day, 16.981 hours/ns, 16.358 timesteps/s
+98.2% CPU use with 1 MPI tasks x no OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
-Pair | 0.16405 | 0.16405 | 0.16405 | 0.0 | 56.61
+Pair | 0.34177 | 0.34177 | 0.34177 | 0.0 | 55.91
Neigh | 0 | 0 | 0 | 0.0 | 0.00
-Comm | 0.00066328 | 0.00066328 | 0.00066328 | 0.0 | 0.23
-Output | 0.037718 | 0.037718 | 0.037718 | 0.0 | 13.02
-Modify | 0.087281 | 0.087281 | 0.087281 | 0.0 | 30.12
-Other | | 7.057e-05 | | | 0.02
+Comm | 0.0013342 | 0.0013342 | 0.0013342 | 0.0 | 0.22
+Output | 0.083583 | 0.083583 | 0.083583 | 0.0 | 13.67
+Modify | 0.18451 | 0.18451 | 0.18451 | 0.0 | 30.18
+Other | | 0.0001087 | | | 0.02
Nlocal: 864 ave 864 max 864 min
Histogram: 1 0 0 0 0 0 0 0 0 0
diff --git a/src/USER-DPD/fix_eos_table_rx.cpp b/src/USER-DPD/fix_eos_table_rx.cpp
index 9a67afd72afa350b8572d0cf0056355df4d5ff20..1dfd58c4658b0372692327562c689f0617951d20 100644
--- a/src/USER-DPD/fix_eos_table_rx.cpp
+++ b/src/USER-DPD/fix_eos_table_rx.cpp
@@ -28,6 +28,12 @@
#define MAXLINE 1024
+#ifdef DBL_EPSILON
+ #define MY_EPSILON (10.0*DBL_EPSILON)
+#else
+ #define MY_EPSILON (10.0*2.220446049250313e-16)
+#endif
+
using namespace LAMMPS_NS;
using namespace FixConst;
@@ -37,17 +43,18 @@ FixEOStableRX::FixEOStableRX(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg), ntables(0), tables(NULL),
tables2(NULL), dHf(NULL), eosSpecies(NULL)
{
- if (narg != 8) error->all(FLERR,"Illegal fix eos/table/rx command");
+ if (narg != 8 && narg != 10) error->all(FLERR,"Illegal fix eos/table/rx command");
restart_peratom = 1;
nevery = 1;
- bool rx_flag = false;
+ rx_flag = false;
+ nspecies = 1;
for (int i = 0; i < modify->nfix; i++)
- if (strncmp(modify->fix[i]->style,"rx",2) == 0) rx_flag = true;
- if (!rx_flag) error->all(FLERR,"FixEOStableRX requires a fix rx command.");
-
- nspecies = atom->nspecies_dpd;
- if(nspecies==0) error->all(FLERR,"There are no rx species specified.");
+ if (strncmp(modify->fix[i]->style,"rx",2) == 0){
+ rx_flag = true;
+ nspecies = atom->nspecies_dpd;
+ if(nspecies==0) error->all(FLERR,"There are no rx species specified.");
+ }
if (strcmp(arg[3],"linear") == 0) tabstyle = LINEAR;
else error->all(FLERR,"Unknown table style in fix eos/table/rx");
@@ -113,8 +120,25 @@ FixEOStableRX::FixEOStableRX(LAMMPS *lmp, int narg, char **arg) :
ntables++;
}
- // Read the Formation Enthalpies
- read_file(arg[7]);
+ // Read the Formation Enthalpies and Correction Coefficients
+ dHf = new double[nspecies];
+ energyCorr = new double[nspecies];
+ tempCorrCoeff = new double[nspecies];
+ moleculeCorrCoeff= new double[nspecies];
+ for (int ii=0; ii<nspecies; ii++){
+ dHf[ii] = 0.0;
+ energyCorr[ii] = 0.0;
+ tempCorrCoeff[ii] = 0.0;
+ moleculeCorrCoeff[ii] = 0.0;
+ }
+
+ if(rx_flag) read_file(arg[7]);
+ else dHf[0] = atof(arg[7]);
+
+ if(narg==10){
+ energyCorr[0] = atof(arg[8]);
+ tempCorrCoeff[0] = atof(arg[9]);
+ }
comm_forward = 3;
comm_reverse = 2;
@@ -136,6 +160,9 @@ FixEOStableRX::~FixEOStableRX()
delete [] dHf;
delete [] eosSpecies;
+ delete [] energyCorr;
+ delete [] tempCorrCoeff;
+ delete [] moleculeCorrCoeff;
}
/* ---------------------------------------------------------------------- */
@@ -269,10 +296,9 @@ void FixEOStableRX::end_of_step()
void FixEOStableRX::read_file(char *file)
{
- int params_per_line = 2;
- char **words = new char*[params_per_line+1];
-
- dHf = new double[nspecies];
+ int min_params_per_line = 2;
+ int max_params_per_line = 5;
+ char **words = new char*[max_params_per_line+1];
// open file on proc 0
@@ -313,7 +339,7 @@ void FixEOStableRX::read_file(char *file)
// concatenate additional lines until have params_per_line words
- while (nwords < params_per_line) {
+ while (nwords < min_params_per_line) {
n = strlen(line);
if (comm->me == 0) {
ptr = fgets(&line[n],MAXLINE-n,fp);
@@ -330,7 +356,7 @@ void FixEOStableRX::read_file(char *file)
nwords = atom->count_words(line);
}
- if (nwords != params_per_line)
+ if (nwords != min_params_per_line && nwords != max_params_per_line)
error->all(FLERR,"Incorrect format in eos table/rx potential file");
// words = ptrs to all words in line
@@ -342,8 +368,14 @@ void FixEOStableRX::read_file(char *file)
for (ispecies = 0; ispecies < nspecies; ispecies++)
if (strcmp(words[0],&atom->dname[ispecies][0]) == 0) break;
- if (ispecies < nspecies)
+ if (ispecies < nspecies){
dHf[ispecies] = atof(words[1]);
+ if(nwords > min_params_per_line+1){
+ energyCorr[ispecies] = atof(words[2]);
+ tempCorrCoeff[ispecies] = atof(words[3]);
+ moleculeCorrCoeff[ispecies] = atof(words[4]);
+ }
+ }
}
delete [] words;
@@ -545,27 +577,33 @@ void FixEOStableRX::param_extract(Table *tb, char *line)
error->one(FLERR,"Invalid keyword in fix eos/table/rx parameters");
word = strtok(NULL," \t\n\r\f");
- while (word) {
- for (ispecies = 0; ispecies < nspecies; ispecies++)
- if (strcmp(word,&atom->dname[ispecies][0]) == 0){
- eosSpecies[ncolumn] = ispecies;
- ncolumn++;
- break;
+ if(rx_flag){
+ while (word) {
+ for (ispecies = 0; ispecies < nspecies; ispecies++)
+ if (strcmp(word,&atom->dname[ispecies][0]) == 0){
+ eosSpecies[ncolumn] = ispecies;
+ ncolumn++;
+ break;
+ }
+ if (ispecies == nspecies){
+ printf("name=%s not found in species list\n",word);
+ error->one(FLERR,"Invalid keyword in fix eos/table/rx parameters");
}
- if (ispecies == nspecies){
- printf("name=%s not found in species list\n",word);
- error->one(FLERR,"Invalid keyword in fix eos/table/rx parameters");
+ word = strtok(NULL," \t\n\r\f");
}
- word = strtok(NULL," \t\n\r\f");
- }
- for (int icolumn = 0; icolumn < ncolumn; icolumn++)
- if(eosSpecies[icolumn]==-1)
- error->one(FLERR,"EOS data is missing from fix eos/table/rx tabe");
- if(ncolumn != nspecies){
- printf("ncolumns=%d nspecies=%d\n",ncolumn,nspecies);
- error->one(FLERR,"The number of columns in fix eos/table/rx does not match the number of species");
+ for (int icolumn = 0; icolumn < ncolumn; icolumn++)
+ if(eosSpecies[icolumn]==-1)
+ error->one(FLERR,"EOS data is missing from fix eos/table/rx tabe");
+ if(ncolumn != nspecies){
+ printf("ncolumns=%d nspecies=%d\n",ncolumn,nspecies);
+ error->one(FLERR,"The number of columns in fix eos/table/rx does not match the number of species");
+ }
+ } else {
+ eosSpecies[0] = 0;
+ ncolumn++;
}
+
if (tb->ninput == 0) error->one(FLERR,"fix eos/table/rx parameters did not set N");
}
@@ -653,11 +691,27 @@ double FixEOStableRX::splint(double *xa, double *ya, double *y2a, int n, double
void FixEOStableRX::energy_lookup(int id, double thetai, double &ui)
{
- int itable;
- double fraction, uTmp, nTotal;
+ int itable, nPG;
+ double fraction, uTmp, nMolecules, nTotal, nTotalPG;
+ double tolerance = 1.0e-10;
ui = 0.0;
nTotal = 0.0;
+ nTotalPG = 0.0;
+ nPG = 0;
+
+ if(rx_flag){
+ for(int ispecies=0;ispecies<nspecies;ispecies++){
+ nTotal += atom->dvector[ispecies][id];
+ if(fabs(moleculeCorrCoeff[ispecies]) > tolerance){
+ nPG++;
+ nTotalPG += atom->dvector[ispecies][id];
+ }
+ }
+ } else {
+ nTotal = 1.0;
+ }
+
for(int ispecies=0;ispecies<nspecies;ispecies++){
Table *tb = &tables[ispecies];
thetai = MAX(thetai,tb->lo);
@@ -669,9 +723,13 @@ void FixEOStableRX::energy_lookup(int id, double thetai, double &ui)
uTmp = tb->e[itable] + fraction*tb->de[itable];
uTmp += dHf[ispecies];
- // mol fraction form:
- ui += atom->dvector[ispecies][id]*uTmp;
- nTotal += atom->dvector[ispecies][id];
+ uTmp += tempCorrCoeff[ispecies]*thetai; // temperature correction
+ uTmp += energyCorr[ispecies]; // energy correction
+ if(nPG > 0) ui += moleculeCorrCoeff[ispecies]*nTotalPG/double(nPG); // molecule correction
+
+ if(rx_flag) nMolecules = atom->dvector[ispecies][id];
+ else nMolecules = 1.0;
+ ui += nMolecules*uTmp;
}
}
ui = ui - double(nTotal+1.5)*force->boltz*thetai;
@@ -690,6 +748,7 @@ void FixEOStableRX::temperature_lookup(int id, double ui, double &thetai)
double maxit = 100;
double temp;
double delta = 0.001;
+ double tolerance = 1.0e-10;
// Store the current thetai in t1
t1 = MAX(thetai,tb->lo);
@@ -713,7 +772,7 @@ void FixEOStableRX::temperature_lookup(int id, double ui, double &thetai)
// Apply the Secant Method
for(it=0; it<maxit; it++){
- if(fabs(f2-f1)<1e-15){
+ if(fabs(f2-f1) < MY_EPSILON){
if(isnan(f1) || isnan(f2)) error->one(FLERR,"NaN detected in secant solver.");
temp = t1;
temp = MAX(temp,tb->lo);
@@ -724,7 +783,7 @@ void FixEOStableRX::temperature_lookup(int id, double ui, double &thetai)
break;
}
temp = t2 - f2*(t2-t1)/(f2-f1);
- if(fabs(temp-t2) < 1e-6) break;
+ if(fabs(temp-t2) < tolerance) break;
f1 = f2;
t1 = t2;
t2 = temp;
diff --git a/src/USER-DPD/fix_eos_table_rx.h b/src/USER-DPD/fix_eos_table_rx.h
index 078cf1e2e1c71ea38eacee620ccc98f4a0123290..8c26d133a52e364bf34bce6c4d7a7b9009adde57 100644
--- a/src/USER-DPD/fix_eos_table_rx.h
+++ b/src/USER-DPD/fix_eos_table_rx.h
@@ -67,7 +67,7 @@ class FixEOStableRX : public Fix {
void read_file(char *);
- double *dHf;
+ double *dHf,*energyCorr,*tempCorrCoeff,*moleculeCorrCoeff;
int pack_reverse_comm(int, int, double *);
void unpack_reverse_comm(int, int *, double *);
@@ -76,6 +76,7 @@ class FixEOStableRX : public Fix {
int *eosSpecies;
int ncolumn;
+ bool rx_flag;
};
}
diff --git a/src/USER-DPD/fix_rx.cpp b/src/USER-DPD/fix_rx.cpp
index e52c032281ede71501459b9986cb514699859d05..0ec814f9a6c22d8cc4ec2b13766697b4cb92cb88 100644
--- a/src/USER-DPD/fix_rx.cpp
+++ b/src/USER-DPD/fix_rx.cpp
@@ -45,6 +45,12 @@ enum{LUCY};
#define MAXLINE 1024
#define DELTA 4
+#ifdef DBL_EPSILON
+ #define MY_EPSILON (10.0*DBL_EPSILON)
+#else
+ #define MY_EPSILON (10.0*2.220446049250313e-16)
+#endif
+
#define SparseKinetics_enableIntegralReactions (true)
#define SparseKinetics_invalidIndex (-1)
@@ -690,7 +696,6 @@ void FixRX::pre_force(int vflag)
int *mask = atom->mask;
double *dpdTheta = atom->dpdTheta;
int newton_pair = force->newton_pair;
- int ii;
double theta;
if(localTempFlag){
@@ -993,9 +998,9 @@ void FixRX::rk4(int id, double *rwork)
// Store the solution back in atom->dvector.
for (int ispecies = 0; ispecies < nspecies; ispecies++){
- if(y[ispecies] < -1.0e-10)
- error->one(FLERR,"Computed concentration in RK4 solver is < -1.0e-10");
- else if(y[ispecies] < 1e-15)
+ if(y[ispecies] < -MY_EPSILON)
+ error->one(FLERR,"Computed concentration in RK4 solver is < -10*DBL_EPSILON");
+ else if(y[ispecies] < MY_EPSILON)
y[ispecies] = 0.0;
atom->dvector[ispecies][id] = y[ispecies];
}
@@ -1512,7 +1517,7 @@ void FixRX::rkf45(int id, double *rwork)
for (int ispecies = 0; ispecies < nspecies; ispecies++){
if(y[ispecies] < -1.0e-10)
error->one(FLERR,"Computed concentration in RKF45 solver is < -1.0e-10");
- else if(y[ispecies] < 1e-20)
+ else if(y[ispecies] < MY_EPSILON)
y[ispecies] = 0.0;
atom->dvector[ispecies][id] = y[ispecies];
}
diff --git a/src/USER-DPD/pair_exp6_rx.cpp b/src/USER-DPD/pair_exp6_rx.cpp
index a1add9a4d2d60ad0a64ca7572fe641bcc66b098f..bf038dd3d702f579cbd7a84dfba2a285b446a248 100644
--- a/src/USER-DPD/pair_exp6_rx.cpp
+++ b/src/USER-DPD/pair_exp6_rx.cpp
@@ -35,6 +35,12 @@ using namespace MathSpecial;
#define MAXLINE 1024
#define DELTA 4
+#ifdef DBL_EPSILON
+ #define MY_EPSILON (10.0*DBL_EPSILON)
+#else
+ #define MY_EPSILON (10.0*2.220446049250313e-16)
+#endif
+
#define oneFluidApproxParameter (-1)
#define isOneFluidApprox(_site) ( (_site) == oneFluidApproxParameter )
@@ -47,17 +53,17 @@ using namespace MathSpecial;
struct PairExp6ParamDataType
{
int n;
- double *epsilon1, *alpha1, *rm1, *fraction1,
- *epsilon2, *alpha2, *rm2, *fraction2,
- *epsilonOld1, *alphaOld1, *rmOld1, *fractionOld1,
- *epsilonOld2, *alphaOld2, *rmOld2, *fractionOld2;
+ double *epsilon1, *alpha1, *rm1, *mixWtSite1,
+ *epsilon2, *alpha2, *rm2, *mixWtSite2,
+ *epsilonOld1, *alphaOld1, *rmOld1, *mixWtSite1old,
+ *epsilonOld2, *alphaOld2, *rmOld2, *mixWtSite2old;
// Default constructor -- nullify everything.
PairExp6ParamDataType(void)
- : n(0), epsilon1(NULL), alpha1(NULL), rm1(NULL), fraction1(NULL),
- epsilon2(NULL), alpha2(NULL), rm2(NULL), fraction2(NULL),
- epsilonOld1(NULL), alphaOld1(NULL), rmOld1(NULL), fractionOld1(NULL),
- epsilonOld2(NULL), alphaOld2(NULL), rmOld2(NULL), fractionOld2(NULL)
+ : n(0), epsilon1(NULL), alpha1(NULL), rm1(NULL), mixWtSite1(NULL),
+ epsilon2(NULL), alpha2(NULL), rm2(NULL), mixWtSite2(NULL),
+ epsilonOld1(NULL), alphaOld1(NULL), rmOld1(NULL), mixWtSite1old(NULL),
+ epsilonOld2(NULL), alphaOld2(NULL), rmOld2(NULL), mixWtSite2old(NULL)
{}
};
@@ -71,6 +77,7 @@ PairExp6rx::PairExp6rx(LAMMPS *lmp) : Pair(lmp)
nparams = maxparam = 0;
params = NULL;
mol2param = NULL;
+ fractionalWeighting = true;
}
/* ---------------------------------------------------------------------- */
@@ -89,6 +96,11 @@ PairExp6rx::~PairExp6rx()
memory->destroy(cutsq);
memory->destroy(cut);
}
+ if(scalingFlag == POLYNOMIAL){
+ memory->destroy(coeffAlpha);
+ memory->destroy(coeffEps);
+ memory->destroy(coeffRm);
+ }
}
/* ---------------------------------------------------------------------- */
@@ -130,10 +142,10 @@ void PairExp6rx::compute(int eflag, int vflag)
double epsilon2_j,alpha2_j,rm2_j;
double evdwlOldEXP6_12, evdwlOldEXP6_21, fpairOldEXP6_12, fpairOldEXP6_21;
double evdwlEXP6_12, evdwlEXP6_21;
- double fractionOld1_i, fractionOld1_j;
- double fractionOld2_i, fractionOld2_j;
- double fraction1_i, fraction1_j;
- double fraction2_i, fraction2_j;
+ double mixWtSite1old_i, mixWtSite1old_j;
+ double mixWtSite2old_i, mixWtSite2old_j;
+ double mixWtSite1_i, mixWtSite1_j;
+ double mixWtSite2_i, mixWtSite2_j;
double *uCG = atom->uCG;
double *uCGnew = atom->uCGnew;
@@ -153,38 +165,38 @@ void PairExp6rx::compute(int eflag, int vflag)
memory->create( PairExp6ParamData.epsilon1 , np_total, "PairExp6ParamData.epsilon1");
memory->create( PairExp6ParamData.alpha1 , np_total, "PairExp6ParamData.alpha1");
memory->create( PairExp6ParamData.rm1 , np_total, "PairExp6ParamData.rm1");
- memory->create( PairExp6ParamData.fraction1 , np_total, "PairExp6ParamData.fraction1");
+ memory->create( PairExp6ParamData.mixWtSite1 , np_total, "PairExp6ParamData.mixWtSite1");
memory->create( PairExp6ParamData.epsilon2 , np_total, "PairExp6ParamData.epsilon2");
memory->create( PairExp6ParamData.alpha2 , np_total, "PairExp6ParamData.alpha2");
memory->create( PairExp6ParamData.rm2 , np_total, "PairExp6ParamData.rm2");
- memory->create( PairExp6ParamData.fraction2 , np_total, "PairExp6ParamData.fraction2");
+ memory->create( PairExp6ParamData.mixWtSite2 , np_total, "PairExp6ParamData.mixWtSite2");
memory->create( PairExp6ParamData.epsilonOld1 , np_total, "PairExp6ParamData.epsilonOld1");
memory->create( PairExp6ParamData.alphaOld1 , np_total, "PairExp6ParamData.alphaOld1");
memory->create( PairExp6ParamData.rmOld1 , np_total, "PairExp6ParamData.rmOld1");
- memory->create( PairExp6ParamData.fractionOld1 , np_total, "PairExp6ParamData.fractionOld1");
+ memory->create( PairExp6ParamData.mixWtSite1old , np_total, "PairExp6ParamData.mixWtSite1old");
memory->create( PairExp6ParamData.epsilonOld2 , np_total, "PairExp6ParamData.epsilonOld2");
memory->create( PairExp6ParamData.alphaOld2 , np_total, "PairExp6ParamData.alphaOld2");
memory->create( PairExp6ParamData.rmOld2 , np_total, "PairExp6ParamData.rmOld2");
- memory->create( PairExp6ParamData.fractionOld2 , np_total, "PairExp6ParamData.fractionOld2");
+ memory->create( PairExp6ParamData.mixWtSite2old , np_total, "PairExp6ParamData.mixWtSite2old");
for (i = 0; i < np_total; ++i)
{
- getParamsEXP6 (i, PairExp6ParamData.epsilon1[i],
+ getMixingWeights (i, PairExp6ParamData.epsilon1[i],
PairExp6ParamData.alpha1[i],
PairExp6ParamData.rm1[i],
- PairExp6ParamData.fraction1[i],
+ PairExp6ParamData.mixWtSite1[i],
PairExp6ParamData.epsilon2[i],
PairExp6ParamData.alpha2[i],
PairExp6ParamData.rm2[i],
- PairExp6ParamData.fraction2[i],
+ PairExp6ParamData.mixWtSite2[i],
PairExp6ParamData.epsilonOld1[i],
PairExp6ParamData.alphaOld1[i],
PairExp6ParamData.rmOld1[i],
- PairExp6ParamData.fractionOld1[i],
+ PairExp6ParamData.mixWtSite1old[i],
PairExp6ParamData.epsilonOld2[i],
PairExp6ParamData.alphaOld2[i],
PairExp6ParamData.rmOld2[i],
- PairExp6ParamData.fractionOld2[i]);
+ PairExp6ParamData.mixWtSite2old[i]);
}
}
@@ -208,19 +220,19 @@ void PairExp6rx::compute(int eflag, int vflag)
epsilon1_i = PairExp6ParamData.epsilon1[i];
alpha1_i = PairExp6ParamData.alpha1[i];
rm1_i = PairExp6ParamData.rm1[i];
- fraction1_i = PairExp6ParamData.fraction1[i];
+ mixWtSite1_i = PairExp6ParamData.mixWtSite1[i];
epsilon2_i = PairExp6ParamData.epsilon2[i];
alpha2_i = PairExp6ParamData.alpha2[i];
rm2_i = PairExp6ParamData.rm2[i];
- fraction2_i = PairExp6ParamData.fraction2[i];
+ mixWtSite2_i = PairExp6ParamData.mixWtSite2[i];
epsilonOld1_i = PairExp6ParamData.epsilonOld1[i];
alphaOld1_i = PairExp6ParamData.alphaOld1[i];
rmOld1_i = PairExp6ParamData.rmOld1[i];
- fractionOld1_i = PairExp6ParamData.fractionOld1[i];
+ mixWtSite1old_i = PairExp6ParamData.mixWtSite1old[i];
epsilonOld2_i = PairExp6ParamData.epsilonOld2[i];
alphaOld2_i = PairExp6ParamData.alphaOld2[i];
rmOld2_i = PairExp6ParamData.rmOld2[i];
- fractionOld2_i = PairExp6ParamData.fractionOld2[i];
+ mixWtSite2old_i = PairExp6ParamData.mixWtSite2old[i];
}
for (jj = 0; jj < jnum; jj++) {
@@ -255,19 +267,19 @@ void PairExp6rx::compute(int eflag, int vflag)
epsilon1_j = PairExp6ParamData.epsilon1[j];
alpha1_j = PairExp6ParamData.alpha1[j];
rm1_j = PairExp6ParamData.rm1[j];
- fraction1_j = PairExp6ParamData.fraction1[j];
+ mixWtSite1_j = PairExp6ParamData.mixWtSite1[j];
epsilon2_j = PairExp6ParamData.epsilon2[j];
alpha2_j = PairExp6ParamData.alpha2[j];
rm2_j = PairExp6ParamData.rm2[j];
- fraction2_j = PairExp6ParamData.fraction2[j];
+ mixWtSite2_j = PairExp6ParamData.mixWtSite2[j];
epsilonOld1_j = PairExp6ParamData.epsilonOld1[j];
alphaOld1_j = PairExp6ParamData.alphaOld1[j];
rmOld1_j = PairExp6ParamData.rmOld1[j];
- fractionOld1_j = PairExp6ParamData.fractionOld1[j];
+ mixWtSite1old_j = PairExp6ParamData.mixWtSite1old[j];
epsilonOld2_j = PairExp6ParamData.epsilonOld2[j];
alphaOld2_j = PairExp6ParamData.alphaOld2[j];
rmOld2_j = PairExp6ParamData.rmOld2[j];
- fractionOld2_j = PairExp6ParamData.fractionOld2[j];
+ mixWtSite2old_j = PairExp6ParamData.mixWtSite2old[j];
}
// A2. Apply Lorentz-Berthelot mixing rules for the i-j pair
@@ -368,9 +380,9 @@ void PairExp6rx::compute(int eflag, int vflag)
}
if (isite1 == isite2)
- evdwlOld = sqrt(fractionOld1_i*fractionOld2_j)*evdwlOldEXP6_12;
+ evdwlOld = sqrt(mixWtSite1old_i*mixWtSite2old_j)*evdwlOldEXP6_12;
else
- evdwlOld = sqrt(fractionOld1_i*fractionOld2_j)*evdwlOldEXP6_12 + sqrt(fractionOld2_i*fractionOld1_j)*evdwlOldEXP6_21;
+ evdwlOld = sqrt(mixWtSite1old_i*mixWtSite2old_j)*evdwlOldEXP6_12 + sqrt(mixWtSite2old_i*mixWtSite1old_j)*evdwlOldEXP6_21;
evdwlOld *= factor_lj;
@@ -451,8 +463,8 @@ void PairExp6rx::compute(int eflag, int vflag)
//
// Apply Mixing Rule to get the overall force for the CG pair
//
- if (isite1 == isite2) fpair = sqrt(fractionOld1_i*fractionOld2_j)*fpairOldEXP6_12;
- else fpair = sqrt(fractionOld1_i*fractionOld2_j)*fpairOldEXP6_12 + sqrt(fractionOld2_i*fractionOld1_j)*fpairOldEXP6_21;
+ if (isite1 == isite2) fpair = sqrt(mixWtSite1old_i*mixWtSite2old_j)*fpairOldEXP6_12;
+ else fpair = sqrt(mixWtSite1old_i*mixWtSite2old_j)*fpairOldEXP6_12 + sqrt(mixWtSite2old_i*mixWtSite1old_j)*fpairOldEXP6_21;
f[i][0] += delx*fpair;
f[i][1] += dely*fpair;
@@ -463,8 +475,8 @@ void PairExp6rx::compute(int eflag, int vflag)
f[j][2] -= delz*fpair;
}
- if (isite1 == isite2) evdwl = sqrt(fraction1_i*fraction2_j)*evdwlEXP6_12;
- else evdwl = sqrt(fraction1_i*fraction2_j)*evdwlEXP6_12 + sqrt(fraction2_i*fraction1_j)*evdwlEXP6_21;
+ if (isite1 == isite2) evdwl = sqrt(mixWtSite1_i*mixWtSite2_j)*evdwlEXP6_12;
+ else evdwl = sqrt(mixWtSite1_i*mixWtSite2_j)*evdwlEXP6_12 + sqrt(mixWtSite2_i*mixWtSite1_j)*evdwlEXP6_21;
evdwl *= factor_lj;
uCGnew[i] += 0.5*evdwl;
@@ -484,19 +496,19 @@ void PairExp6rx::compute(int eflag, int vflag)
if (PairExp6ParamData.epsilon1 ) memory->destroy(PairExp6ParamData.epsilon1);
if (PairExp6ParamData.alpha1 ) memory->destroy(PairExp6ParamData.alpha1);
if (PairExp6ParamData.rm1 ) memory->destroy(PairExp6ParamData.rm1);
- if (PairExp6ParamData.fraction1 ) memory->destroy(PairExp6ParamData.fraction1);
+ if (PairExp6ParamData.mixWtSite1 ) memory->destroy(PairExp6ParamData.mixWtSite1);
if (PairExp6ParamData.epsilon2 ) memory->destroy(PairExp6ParamData.epsilon2);
if (PairExp6ParamData.alpha2 ) memory->destroy(PairExp6ParamData.alpha2);
if (PairExp6ParamData.rm2 ) memory->destroy(PairExp6ParamData.rm2);
- if (PairExp6ParamData.fraction2 ) memory->destroy(PairExp6ParamData.fraction2);
+ if (PairExp6ParamData.mixWtSite2 ) memory->destroy(PairExp6ParamData.mixWtSite2);
if (PairExp6ParamData.epsilonOld1 ) memory->destroy(PairExp6ParamData.epsilonOld1);
if (PairExp6ParamData.alphaOld1 ) memory->destroy(PairExp6ParamData.alphaOld1);
if (PairExp6ParamData.rmOld1 ) memory->destroy(PairExp6ParamData.rmOld1);
- if (PairExp6ParamData.fractionOld1) memory->destroy(PairExp6ParamData.fractionOld1);
+ if (PairExp6ParamData.mixWtSite1old) memory->destroy(PairExp6ParamData.mixWtSite1old);
if (PairExp6ParamData.epsilonOld2 ) memory->destroy(PairExp6ParamData.epsilonOld2);
if (PairExp6ParamData.alphaOld2 ) memory->destroy(PairExp6ParamData.alphaOld2);
if (PairExp6ParamData.rmOld2 ) memory->destroy(PairExp6ParamData.rmOld2);
- if (PairExp6ParamData.fractionOld2) memory->destroy(PairExp6ParamData.fractionOld2);
+ if (PairExp6ParamData.mixWtSite2old) memory->destroy(PairExp6ParamData.mixWtSite2old);
}
}
@@ -526,10 +538,20 @@ void PairExp6rx::allocate()
void PairExp6rx::settings(int narg, char **arg)
{
- if (narg != 1) error->all(FLERR,"Illegal pair_style command");
+ if (narg < 1) error->all(FLERR,"Illegal pair_style command");
cut_global = force->numeric(FLERR,arg[0]);
+ // optional keywords
+
+ int iarg = 1;
+ while (iarg < narg) {
+ if (strcmp(arg[iarg],"fractional") == 0) fractionalWeighting = true;
+ else if (strcmp(arg[iarg],"molecular") == 0) fractionalWeighting = false;
+ else error->all(FLERR,"Illegal pair_style command");
+ iarg++;
+ }
+
if (allocated) {
int i,j;
for (i = 1; i <= atom->ntypes; i++)
@@ -547,7 +569,7 @@ void PairExp6rx::settings(int narg, char **arg)
void PairExp6rx::coeff(int narg, char **arg)
{
- if (narg < 7 || narg > 8) error->all(FLERR,"Incorrect args for pair coefficients");
+ if (narg < 6 || narg > 9) error->all(FLERR,"Incorrect args for pair coefficients");
bool rx_flag = false;
for (int i = 0; i < modify->nfix; i++)
@@ -624,21 +646,36 @@ void PairExp6rx::coeff(int narg, char **arg)
params[iparam].potentialType = exp6PotentialType;
else
error->all(FLERR,"params[].potential type unknown");
-
- //printf("params[%d].name= %s ispecies= %d potential= %s potentialType= %d\n", iparam, params[iparam].name, params[iparam].ispecies, params[iparam].potential, params[iparam].potentialType);
}
}
delete[] site1;
delete[] site2;
site1 = site2 = NULL;
- fuchslinR = force->numeric(FLERR,arg[5]);
- fuchslinEpsilon = force->numeric(FLERR,arg[6]);
-
setup();
double cut_one = cut_global;
- if (narg == 8) cut_one = force->numeric(FLERR,arg[7]);
+ if (strcmp(arg[5],"exponent") == 0){
+ scalingFlag = EXPONENT;
+ exponentR = force->numeric(FLERR,arg[6]);
+ exponentEpsilon = force->numeric(FLERR,arg[7]);
+ if (narg > 9) error->all(FLERR,"Incorrect args for pair coefficients");
+ if (narg == 9) cut_one = force->numeric(FLERR,arg[8]);
+ } else if (strcmp(arg[5],"polynomial") == 0){
+ scalingFlag = POLYNOMIAL;
+ memory->create(coeffAlpha,6,"pair:coeffAlpha");
+ memory->create(coeffEps,6,"pair:coeffEps");
+ memory->create(coeffRm,6,"pair:coeffRm");
+ read_file2(arg[6]);
+ if (narg > 8) error->all(FLERR,"Incorrect args for pair coefficients");
+ if (narg == 8) cut_one = force->numeric(FLERR,arg[7]);
+ } else if (strcmp(arg[5],"none") == 0){
+ scalingFlag = NONE;
+ if (narg > 7) error->all(FLERR,"Incorrect args for pair coefficients");
+ if (narg == 7) cut_one = force->numeric(FLERR,arg[6]);
+ } else {
+ error->all(FLERR,"Incorrect args for pair coefficients");
+ }
int count = 0;
for (int i = ilo; i <= ihi; i++) {
@@ -780,6 +817,95 @@ void PairExp6rx::read_file(char *file)
/* ---------------------------------------------------------------------- */
+void PairExp6rx::read_file2(char *file)
+{
+ int params_per_line = 7;
+ char **words = new char*[params_per_line+1];
+
+ // open file on proc 0
+
+ FILE *fp;
+ fp = NULL;
+ if (comm->me == 0) {
+ fp = fopen(file,"r");
+ if (fp == NULL) {
+ char str[128];
+ sprintf(str,"Cannot open polynomial file %s",file);
+ error->one(FLERR,str);
+ }
+ }
+
+ // one set of params can span multiple lines
+ int n,nwords;
+ 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 polynomial 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;
+
+ if (strcmp(words[0],"alpha") == 0){
+ for (int ii=1; ii<params_per_line; ii++)
+ coeffAlpha[ii-1] = atof(words[ii]);
+ }
+ if (strcmp(words[0],"epsilon") == 0){
+ for (int ii=1; ii<params_per_line; ii++)
+ coeffEps[ii-1] = atof(words[ii]);
+ }
+ if (strcmp(words[0],"rm") == 0){
+ for (int ii=1; ii<params_per_line; ii++)
+ coeffRm[ii-1] = atof(words[ii]);
+ }
+ }
+
+ delete [] words;
+}
+
+/* ---------------------------------------------------------------------- */
+
void PairExp6rx::setup()
{
int i,j,n;
@@ -877,7 +1003,7 @@ void PairExp6rx::read_restart_settings(FILE *fp)
/* ---------------------------------------------------------------------- */
-void PairExp6rx::getParamsEXP6(int id,double &epsilon1,double &alpha1,double &rm1, double &fraction1,double &epsilon2,double &alpha2,double &rm2,double &fraction2,double &epsilon1_old,double &alpha1_old,double &rm1_old, double &fraction1_old,double &epsilon2_old,double &alpha2_old,double &rm2_old,double &fraction2_old) const
+void PairExp6rx::getMixingWeights(int id,double &epsilon1,double &alpha1,double &rm1, double &mixWtSite1,double &epsilon2,double &alpha2,double &rm2,double &mixWtSite2,double &epsilon1_old,double &alpha1_old,double &rm1_old, double &mixWtSite1old,double &epsilon2_old,double &alpha2_old,double &rm2_old,double &mixWtSite2old) const
{
int iparam, jparam;
double rmi, rmj, rmij, rm3ij;
@@ -885,11 +1011,16 @@ void PairExp6rx::getParamsEXP6(int id,double &epsilon1,double &alpha1,double &rm
double alphai, alphaj, alphaij;
double epsilon_old, rm3_old, alpha_old;
double epsilon, rm3, alpha;
- double fractionOFA, fractionOFA_old;
- double nTotalOFA, nTotalOFA_old;
- double nTotal, nTotal_old;
double xMolei, xMolej, xMolei_old, xMolej_old;
+ double fractionOFAold, fractionOFA;
+ double fractionOld1, fraction1;
+ double fractionOld2, fraction2;
+ double nMoleculesOFAold, nMoleculesOFA;
+ double nMoleculesOld1, nMolecules1;
+ double nMoleculesOld2, nMolecules2;
+ double nTotal, nTotalOld;
+
rm3 = 0.0;
epsilon = 0.0;
alpha = 0.0;
@@ -897,32 +1028,32 @@ void PairExp6rx::getParamsEXP6(int id,double &epsilon1,double &alpha1,double &rm
rm3_old = 0.0;
alpha_old = 0.0;
fractionOFA = 0.0;
- fractionOFA_old = 0.0;
- nTotalOFA = 0.0;
- nTotalOFA_old = 0.0;
+ fractionOFAold = 0.0;
+ nMoleculesOFA = 0.0;
+ nMoleculesOFAold = 0.0;
nTotal = 0.0;
- nTotal_old = 0.0;
+ nTotalOld = 0.0;
// Compute the total number of molecules in the old and new CG particle as well as the total number of molecules in the fluid portion of the old and new CG particle
for (int ispecies = 0; ispecies < nspecies; ispecies++){
nTotal += atom->dvector[ispecies][id];
- nTotal_old += atom->dvector[ispecies+nspecies][id];
+ nTotalOld += atom->dvector[ispecies+nspecies][id];
iparam = mol2param[ispecies];
if (iparam < 0 || params[iparam].potentialType != exp6PotentialType ) continue;
if (isOneFluidApprox(isite1) || isOneFluidApprox(isite2)) {
if (isite1 == params[iparam].ispecies || isite2 == params[iparam].ispecies) continue;
- nTotalOFA_old += atom->dvector[ispecies+nspecies][id];
- nTotalOFA += atom->dvector[ispecies][id];
+ nMoleculesOFAold += atom->dvector[ispecies+nspecies][id];
+ nMoleculesOFA += atom->dvector[ispecies][id];
}
}
- if(nTotal < 1e-8 || nTotal_old < 1e-8)
- error->all(FLERR,"The number of molecules in CG particle is less than 1e-8.");
+ if(nTotal < MY_EPSILON || nTotalOld < MY_EPSILON)
+ error->all(FLERR,"The number of molecules in CG particle is less than 10*DBL_EPSILON.");
// Compute the mole fraction of molecules within the fluid portion of the particle (One Fluid Approximation)
- fractionOFA_old = nTotalOFA_old / nTotal_old;
- fractionOFA = nTotalOFA / nTotal;
+ fractionOFAold = nMoleculesOFAold / nTotalOld;
+ fractionOFA = nMoleculesOFA / nTotal;
for (int ispecies = 0; ispecies < nspecies; ispecies++) {
iparam = mol2param[ispecies];
@@ -938,8 +1069,10 @@ void PairExp6rx::getParamsEXP6(int id,double &epsilon1,double &alpha1,double &rm
alpha1 = params[iparam].alpha;
// Compute the mole fraction of Site1
- fraction1_old = atom->dvector[ispecies+nspecies][id]/nTotal_old;
- fraction1 = atom->dvector[ispecies][id]/nTotal;
+ nMoleculesOld1 = atom->dvector[ispecies+nspecies][id];
+ nMolecules1 = atom->dvector[ispecies][id];
+ fractionOld1 = nMoleculesOld1/nTotalOld;
+ fraction1 = nMolecules1/nTotal;
}
// If Site2 matches a pure species, then grab the parameters
@@ -952,7 +1085,9 @@ void PairExp6rx::getParamsEXP6(int id,double &epsilon1,double &alpha1,double &rm
alpha2 = params[iparam].alpha;
// Compute the mole fraction of Site2
- fraction2_old = atom->dvector[ispecies+nspecies][id]/nTotal_old;
+ nMoleculesOld2 = atom->dvector[ispecies+nspecies][id];
+ nMolecules2 = atom->dvector[ispecies][id];
+ fractionOld2 = atom->dvector[ispecies+nspecies][id]/nTotalOld;
fraction2 = atom->dvector[ispecies][id]/nTotal;
}
@@ -962,8 +1097,10 @@ void PairExp6rx::getParamsEXP6(int id,double &epsilon1,double &alpha1,double &rm
rmi = params[iparam].rm;
epsiloni = params[iparam].epsilon;
alphai = params[iparam].alpha;
- xMolei = atom->dvector[ispecies][id]/nTotalOFA;
- xMolei_old = atom->dvector[ispecies+nspecies][id]/nTotalOFA_old;
+ if(nMoleculesOFA<MY_EPSILON) xMolei = 0.0;
+ else xMolei = atom->dvector[ispecies][id]/nMoleculesOFA;
+ if(nMoleculesOFAold<MY_EPSILON) xMolei_old = 0.0;
+ else xMolei_old = atom->dvector[ispecies+nspecies][id]/nMoleculesOFAold;
for (int jspecies = 0; jspecies < nspecies; jspecies++) {
jparam = mol2param[jspecies];
@@ -972,15 +1109,17 @@ void PairExp6rx::getParamsEXP6(int id,double &epsilon1,double &alpha1,double &rm
rmj = params[jparam].rm;
epsilonj = params[jparam].epsilon;
alphaj = params[jparam].alpha;
- xMolej = atom->dvector[jspecies][id]/nTotalOFA;
- xMolej_old = atom->dvector[jspecies+nspecies][id]/nTotalOFA_old;
+ if(nMoleculesOFA<MY_EPSILON) xMolej = 0.0;
+ else xMolej = atom->dvector[jspecies][id]/nMoleculesOFA;
+ if(nMoleculesOFAold<MY_EPSILON) xMolej_old = 0.0;
+ else xMolej_old = atom->dvector[jspecies+nspecies][id]/nMoleculesOFAold;
rmij = (rmi+rmj)/2.0;
rm3ij = rmij*rmij*rmij;
epsilonij = sqrt(epsiloni*epsilonj);
alphaij = sqrt(alphai*alphaj);
- if(fractionOFA_old > 0.0){
+ if(fractionOFAold > 0.0){
rm3_old += xMolei_old*xMolej_old*rm3ij;
epsilon_old += xMolei_old*xMolej_old*rm3ij*epsilonij;
alpha_old += xMolei_old*xMolej_old*rm3ij*epsilonij*alphaij;
@@ -996,7 +1135,7 @@ void PairExp6rx::getParamsEXP6(int id,double &epsilon1,double &alpha1,double &rm
if (isOneFluidApprox(isite1)){
rm1 = cbrt(rm3);
- if(rm1 < 1e-16) {
+ if(rm1 < MY_EPSILON) {
rm1 = 0.0;
epsilon1 = 0.0;
alpha1 = 0.0;
@@ -1004,11 +1143,11 @@ void PairExp6rx::getParamsEXP6(int id,double &epsilon1,double &alpha1,double &rm
epsilon1 = epsilon / rm3;
alpha1 = alpha / epsilon1 / rm3;
}
-
+ nMolecules1 = 1.0-(nTotal-nMoleculesOFA);
fraction1 = fractionOFA;
rm1_old = cbrt(rm3_old);
- if(rm1_old < 1e-16) {
+ if(rm1_old < MY_EPSILON) {
rm1_old = 0.0;
epsilon1_old = 0.0;
alpha1_old = 0.0;
@@ -1016,42 +1155,21 @@ void PairExp6rx::getParamsEXP6(int id,double &epsilon1,double &alpha1,double &rm
epsilon1_old = epsilon_old / rm3_old;
alpha1_old = alpha_old / epsilon1_old / rm3_old;
}
- fraction1_old = fractionOFA_old;
-
- // Fuchslin-Like Exp-6 Scaling
- double powfuch = 0.0;
- if(fuchslinEpsilon < 0.0){
- powfuch = pow(nTotalOFA,-fuchslinEpsilon);
- if(powfuch<1e-15) epsilon1 = 0.0;
- else epsilon1 *= 1.0/powfuch;
-
- powfuch = pow(nTotalOFA_old,-fuchslinEpsilon);
- if(powfuch<1e-15) epsilon1_old = 0.0;
- else epsilon1_old *= 1.0/powfuch;
-
- } else {
- epsilon1 *= pow(nTotalOFA,fuchslinEpsilon);
- epsilon1_old *= pow(nTotalOFA_old,fuchslinEpsilon);
- }
-
- if(fuchslinR < 0.0){
- powfuch = pow(nTotalOFA,-fuchslinR);
- if(powfuch<1e-15) rm1 = 0.0;
- else rm1 *= 1.0/powfuch;
-
- powfuch = pow(nTotalOFA_old,-fuchslinR);
- if(powfuch<1e-15) rm1_old = 0.0;
- else rm1_old *= 1.0/powfuch;
-
- } else {
- rm1 *= pow(nTotalOFA,fuchslinR);
- rm1_old *= pow(nTotalOFA_old,fuchslinR);
+ nMoleculesOld1 = 1.0-(nTotalOld-nMoleculesOFAold);
+ fractionOld1 = fractionOFAold;
+
+ if(scalingFlag == EXPONENT){
+ exponentScaling(nMoleculesOFA,epsilon1,rm1);
+ exponentScaling(nMoleculesOFAold,epsilon1_old,rm1_old);
+ } else if(scalingFlag == POLYNOMIAL){
+ polynomialScaling(nMoleculesOFA,alpha1,epsilon1,rm1);
+ polynomialScaling(nMoleculesOFAold,alpha1_old,epsilon1_old,rm1_old);
}
}
if (isOneFluidApprox(isite2)){
rm2 = cbrt(rm3);
- if(rm2 < 1e-16) {
+ if(rm2 < MY_EPSILON) {
rm2 = 0.0;
epsilon2 = 0.0;
alpha2 = 0.0;
@@ -1059,10 +1177,11 @@ void PairExp6rx::getParamsEXP6(int id,double &epsilon1,double &alpha1,double &rm
epsilon2 = epsilon / rm3;
alpha2 = alpha / epsilon2 / rm3;
}
+ nMolecules2 = 1.0-(nTotal-nMoleculesOFA);
fraction2 = fractionOFA;
rm2_old = cbrt(rm3_old);
- if(rm2_old < 1e-16) {
+ if(rm2_old < MY_EPSILON) {
rm2_old = 0.0;
epsilon2_old = 0.0;
alpha2_old = 0.0;
@@ -1070,64 +1189,96 @@ void PairExp6rx::getParamsEXP6(int id,double &epsilon1,double &alpha1,double &rm
epsilon2_old = epsilon_old / rm3_old;
alpha2_old = alpha_old / epsilon2_old / rm3_old;
}
- fraction2_old = fractionOFA_old;
-
- // Fuchslin-Like Exp-6 Scaling
- double powfuch = 0.0;
- if(fuchslinEpsilon < 0.0){
- powfuch = pow(nTotalOFA,-fuchslinEpsilon);
- if(powfuch<1e-15) epsilon2 = 0.0;
- else epsilon2 *= 1.0/powfuch;
-
- powfuch = pow(nTotalOFA_old,-fuchslinEpsilon);
- if(powfuch<1e-15) epsilon2_old = 0.0;
- else epsilon2_old *= 1.0/powfuch;
-
- } else {
- epsilon2 *= pow(nTotalOFA,fuchslinEpsilon);
- epsilon2_old *= pow(nTotalOFA_old,fuchslinEpsilon);
- }
-
- if(fuchslinR < 0.0){
- powfuch = pow(nTotalOFA,-fuchslinR);
- if(powfuch<1e-15) rm2 = 0.0;
- else rm2 *= 1.0/powfuch;
-
- powfuch = pow(nTotalOFA_old,-fuchslinR);
- if(powfuch<1e-15) rm2_old = 0.0;
- else rm2_old *= 1.0/powfuch;
-
- } else {
- rm2 *= pow(nTotalOFA,fuchslinR);
- rm2_old *= pow(nTotalOFA_old,fuchslinR);
+ nMoleculesOld2 = 1.0-(nTotalOld-nMoleculesOFAold);
+ fractionOld2 = fractionOFAold;
+
+ if(scalingFlag == EXPONENT){
+ exponentScaling(nMoleculesOFA,epsilon2,rm2);
+ exponentScaling(nMoleculesOFAold,epsilon2_old,rm2_old);
+ } else if(scalingFlag == POLYNOMIAL){
+ polynomialScaling(nMoleculesOFA,alpha2,epsilon2,rm2);
+ polynomialScaling(nMoleculesOFAold,alpha2_old,epsilon2_old,rm2_old);
}
}
// Check that no fractions are less than zero
- if(fraction1 < 0.0){
- if(fraction1 < -1.0e-10){
- error->all(FLERR,"Computed fraction less than -1.0e-10");
+ if(fraction1 < 0.0 || nMolecules1 < 0.0){
+ if(fraction1 < -MY_EPSILON || nMolecules1 < -MY_EPSILON){
+ error->all(FLERR,"Computed fraction less than -10*DBL_EPSILON");
}
+ nMolecules1 = 0.0;
fraction1 = 0.0;
}
- if(fraction2 < 0.0){
- if(fraction2 < -1.0e-10){
- error->all(FLERR,"Computed fraction less than -1.0e-10");
+ if(fraction2 < 0.0 || nMolecules2 < 0.0){
+ if(fraction2 < -MY_EPSILON || nMolecules2 < -MY_EPSILON){
+ error->all(FLERR,"Computed fraction less than -10*DBL_EPSILON");
}
+ nMolecules2 = 0.0;
fraction2 = 0.0;
}
- if(fraction1_old < 0.0){
- if(fraction1_old < -1.0e-10){
- error->all(FLERR,"Computed fraction less than -1.0e-10");
+ if(fractionOld1 < 0.0 || nMoleculesOld1 < 0.0){
+ if(fractionOld1 < -MY_EPSILON || nMoleculesOld1 < -MY_EPSILON){
+ error->all(FLERR,"Computed fraction less than -10*DBL_EPSILON");
}
- fraction1_old = 0.0;
+ nMoleculesOld1 = 0.0;
+ fractionOld1 = 0.0;
}
- if(fraction2_old < 0.0){
- if(fraction2_old < -1.0e-10){
- error->all(FLERR,"Computed fraction less than -1.0e-10");
+ if(fractionOld2 < 0.0 || nMoleculesOld2 < 0.0){
+ if(fractionOld2 < -MY_EPSILON || nMoleculesOld2 < -MY_EPSILON){
+ error->all(FLERR,"Computed fraction less than -10*DBL_EPSILON");
}
- fraction2_old = 0.0;
+ nMoleculesOld2 = 0.0;
+ fractionOld2 = 0.0;
}
+
+ if(fractionalWeighting){
+ mixWtSite1old = fractionOld1;
+ mixWtSite1 = fraction1;
+ mixWtSite2old = fractionOld2;
+ mixWtSite2 = fraction2;
+ } else {
+ mixWtSite1old = nMoleculesOld1;
+ mixWtSite1 = nMolecules1;
+ mixWtSite2old = nMoleculesOld2;
+ mixWtSite2 = nMolecules2;
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairExp6rx::exponentScaling(double phi, double &epsilon, double &rm) const
+{
+ double powfuch;
+
+ if(exponentEpsilon < 0.0){
+ powfuch = pow(phi,-exponentEpsilon);
+ if(powfuch<MY_EPSILON) epsilon = 0.0;
+ else epsilon *= 1.0/powfuch;
+ } else {
+ epsilon *= pow(phi,exponentEpsilon);
+ }
+
+ if(exponentR < 0.0){
+ powfuch = pow(phi,-exponentR);
+ if(powfuch<MY_EPSILON) rm = 0.0;
+ else rm *= 1.0/powfuch;
+ } else {
+ rm *= pow(phi,exponentR);
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairExp6rx::polynomialScaling(double phi, double &alpha, double &epsilon, double &rm) const
+{
+ double phi2 = phi*phi;
+ double phi3 = phi2*phi;
+ double phi4 = phi2*phi2;
+ double phi5 = phi2*phi3;
+
+ alpha = (coeffAlpha[0]*phi5 + coeffAlpha[1]*phi4 + coeffAlpha[2]*phi3 + coeffAlpha[3]*phi2 + coeffAlpha[4]*phi + coeffAlpha[5]);
+ epsilon *= (coeffEps[0]*phi5 + coeffEps[1]*phi4 + coeffEps[2]*phi3 + coeffEps[3]*phi2 + coeffEps[4]*phi + coeffEps[5]);
+ rm *= (coeffRm[0]*phi5 + coeffRm[1]*phi4 + coeffRm[2]*phi3 + coeffRm[3]*phi2 + coeffRm[4]*phi + coeffRm[5]);
}
/* ---------------------------------------------------------------------- */
diff --git a/src/USER-DPD/pair_exp6_rx.h b/src/USER-DPD/pair_exp6_rx.h
index dd9fa22a48096f5b8ce3a7b3d07aa3b2b83aa2f2..c16875c96021702ecb532aed737b688daf8d7303 100644
--- a/src/USER-DPD/pair_exp6_rx.h
+++ b/src/USER-DPD/pair_exp6_rx.h
@@ -39,6 +39,7 @@ class PairExp6rx : public Pair {
protected:
enum{LINEAR};
+ enum{NONE,EXPONENT,POLYNOMIAL};
double cut_global;
double **cut;
double **epsilon,**rm,**alpha;
@@ -59,12 +60,18 @@ class PairExp6rx : public Pair {
int nspecies;
void read_file(char *);
+ void read_file2(char *);
void setup();
int isite1, isite2;
char *site1, *site2;
- double fuchslinR, fuchslinEpsilon;
- void getParamsEXP6(int, double &, double &, double &, double &, double &, double &, double &, double &, double &, double &, double &, double &, double &, double &, double &, double &) const;
+ void getMixingWeights(int, double &, double &, double &, double &, double &, double &, double &, double &, double &, double &, double &, double &, double &, double &, double &, double &) const;
+ double exponentR, exponentEpsilon;
+ int scalingFlag;
+ void exponentScaling(double, double &, double &) const;
+ void polynomialScaling(double, double &, double &, double &) const;
+ double *coeffAlpha, *coeffEps, *coeffRm;
+ bool fractionalWeighting;
inline double func_rin(const double &) const;
inline double expValue(const double) const;
@@ -128,7 +135,7 @@ E: Potential file has duplicate entry.
Self-explanatory
-E: The number of molecules in CG particle is less than 1e-8.
+E: The number of molecules in CG particle is less than 10*DBL_EPSILON.
Self-explanatory. Check the species concentrations have been properly set
and check the reaction kinetic solver parameters in fix rx to more for
diff --git a/src/USER-DPD/pair_multi_lucy_rx.cpp b/src/USER-DPD/pair_multi_lucy_rx.cpp
index cd107f1519e8f5918eb6ecc71ad50f4086c207d0..6e20f206d5b81a114e10e3997044bb6b60adb876 100644
--- a/src/USER-DPD/pair_multi_lucy_rx.cpp
+++ b/src/USER-DPD/pair_multi_lucy_rx.cpp
@@ -43,6 +43,12 @@ enum{NONE,RLINEAR,RSQ};
#define MAXLINE 1024
+#ifdef DBL_EPSILON
+ #define MY_EPSILON (10.0*DBL_EPSILON)
+#else
+ #define MY_EPSILON (10.0*2.220446049250313e-16)
+#endif
+
#define oneFluidParameter (-1)
#define isOneFluid(_site) ( (_site) == oneFluidParameter )
@@ -72,6 +78,7 @@ PairMultiLucyRX::PairMultiLucyRX(LAMMPS *lmp) : Pair(lmp),
comm_forward = 1;
comm_reverse = 1;
+ fractionalWeighting = true;
}
/* ---------------------------------------------------------------------- */
@@ -112,9 +119,9 @@ void PairMultiLucyRX::compute(int eflag, int vflag)
int nghost = atom->nghost;
int newton_pair = force->newton_pair;
- double fractionOld1_i,fractionOld1_j;
- double fractionOld2_i,fractionOld2_j;
- double fraction1_i;
+ double mixWtSite1old_i,mixWtSite1old_j;
+ double mixWtSite2old_i,mixWtSite2old_j;
+ double mixWtSite1_i;
double *uCG = atom->uCG;
double *uCGnew = atom->uCGnew;
@@ -124,20 +131,20 @@ void PairMultiLucyRX::compute(int eflag, int vflag)
int jtable;
double *rho = atom->rho;
- double *fractionOld1 = NULL;
- double *fractionOld2 = NULL;
- double *fraction1 = NULL;
- double *fraction2 = NULL;
+ double *mixWtSite1old = NULL;
+ double *mixWtSite2old = NULL;
+ double *mixWtSite1 = NULL;
+ double *mixWtSite2 = NULL;
{
const int ntotal = nlocal + nghost;
- memory->create(fractionOld1, ntotal, "PairMultiLucyRX::fractionOld1");
- memory->create(fractionOld2, ntotal, "PairMultiLucyRX::fractionOld2");
- memory->create(fraction1, ntotal, "PairMultiLucyRX::fraction1");
- memory->create(fraction2, ntotal, "PairMultiLucyRX::fraction2");
+ memory->create(mixWtSite1old, ntotal, "PairMultiLucyRX::mixWtSite1old");
+ memory->create(mixWtSite2old, ntotal, "PairMultiLucyRX::mixWtSite2old");
+ memory->create(mixWtSite1, ntotal, "PairMultiLucyRX::mixWtSite1");
+ memory->create(mixWtSite2, ntotal, "PairMultiLucyRX::mixWtSite2");
for (int i = 0; i < ntotal; ++i)
- getParams(i, fractionOld1[i], fractionOld2[i], fraction1[i], fraction2[i]);
+ getMixingWeights(i, mixWtSite1old[i], mixWtSite2old[i], mixWtSite1[i], mixWtSite2[i]);
}
inum = list->inum;
@@ -162,9 +169,9 @@ void PairMultiLucyRX::compute(int eflag, int vflag)
double fy_i = 0.0;
double fz_i = 0.0;
- fractionOld1_i = fractionOld1[i];
- fractionOld2_i = fractionOld2[i];
- fraction1_i = fraction1[i];
+ mixWtSite1old_i = mixWtSite1old[i];
+ mixWtSite2old_i = mixWtSite2old[i];
+ mixWtSite1_i = mixWtSite1[i];
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
@@ -179,8 +186,8 @@ void PairMultiLucyRX::compute(int eflag, int vflag)
if (rsq < cutsq[itype][jtype]) {
fpair = 0.0;
- fractionOld1_j = fractionOld1[j];
- fractionOld2_j = fractionOld2[j];
+ mixWtSite1old_j = mixWtSite1old[j];
+ mixWtSite2old_j = mixWtSite2old[j];
tb = &tables[tabindex[itype][jtype]];
if (rho[i]*rho[i] < tb->innersq || rho[j]*rho[j] < tb->innersq){
@@ -235,8 +242,8 @@ void PairMultiLucyRX::compute(int eflag, int vflag)
} else error->one(FLERR,"Only LOOKUP and LINEAR table styles have been implemented for pair multi/lucy/rx");
- if (isite1 == isite2) fpair = sqrt(fractionOld1_i*fractionOld2_j)*fpair;
- else fpair = (sqrt(fractionOld1_i*fractionOld2_j) + sqrt(fractionOld2_i*fractionOld1_j))*fpair;
+ if (isite1 == isite2) fpair = sqrt(mixWtSite1old_i*mixWtSite2old_j)*fpair;
+ else fpair = (sqrt(mixWtSite1old_i*mixWtSite2old_j) + sqrt(mixWtSite2old_i*mixWtSite1old_j))*fpair;
fx_i += delx*fpair;
fy_i += dely*fpair;
@@ -268,8 +275,8 @@ void PairMultiLucyRX::compute(int eflag, int vflag)
} else error->one(FLERR,"Only LOOKUP and LINEAR table styles have been implemented for pair multi/lucy/rx");
evdwl *=(pi*cutsq[itype][itype]*cutsq[itype][itype])/84.0;
- evdwlOld = fractionOld1_i*evdwl;
- evdwl = fraction1_i*evdwl;
+ evdwlOld = mixWtSite1old_i*evdwl;
+ evdwl = mixWtSite1_i*evdwl;
uCG[i] += evdwlOld;
uCGnew[i] += evdwl;
@@ -281,10 +288,10 @@ void PairMultiLucyRX::compute(int eflag, int vflag)
if (vflag_fdotr) virial_fdotr_compute();
- memory->destroy(fractionOld1);
- memory->destroy(fractionOld2);
- memory->destroy(fraction1);
- memory->destroy(fraction2);
+ memory->destroy(mixWtSite1old);
+ memory->destroy(mixWtSite2old);
+ memory->destroy(mixWtSite1);
+ memory->destroy(mixWtSite2);
}
/* ----------------------------------------------------------------------
@@ -311,7 +318,7 @@ void PairMultiLucyRX::allocate()
void PairMultiLucyRX::settings(int narg, char **arg)
{
- if (narg != 2) error->all(FLERR,"Illegal pair_style command");
+ if (narg < 2) error->all(FLERR,"Illegal pair_style command");
// new settings
@@ -322,6 +329,16 @@ void PairMultiLucyRX::settings(int narg, char **arg)
tablength = force->inumeric(FLERR,arg[1]);
if (tablength < 2) error->all(FLERR,"Illegal number of pair table entries");
+ // optional keywords
+
+ int iarg = 2;
+ while (iarg < narg) {
+ if (strcmp(arg[iarg],"fractional") == 0) fractionalWeighting = true;
+ else if (strcmp(arg[iarg],"molecular") == 0) fractionalWeighting = false;
+ else error->all(FLERR,"Illegal pair_style command");
+ iarg++;
+ }
+
// delete old tables, since cannot just change settings
for (int m = 0; m < ntables; m++) free_table(&tables[m]);
@@ -928,9 +945,14 @@ void PairMultiLucyRX::computeLocalDensity()
/* ---------------------------------------------------------------------- */
-void PairMultiLucyRX::getParams(int id, double &fractionOld1, double &fractionOld2, double &fraction1, double &fraction2)
+void PairMultiLucyRX::getMixingWeights(int id, double &mixWtSite1old, double &mixWtSite2old, double &mixWtSite1, double &mixWtSite2)
{
- double fractionOld, fraction;
+ double fractionOFAold, fractionOFA;
+ double fractionOld1, fraction1;
+ double fractionOld2, fraction2;
+ double nMoleculesOFAold, nMoleculesOFA;
+ double nMoleculesOld1, nMolecules1;
+ double nMoleculesOld2, nMolecules2;
double nTotal, nTotalOld;
nTotal = 0.0;
@@ -941,32 +963,56 @@ void PairMultiLucyRX::getParams(int id, double &fractionOld1, double &fractionOl
}
if (isOneFluid(isite1) == false){
- fractionOld1 = atom->dvector[isite1+nspecies][id]/nTotalOld;
- fraction1 = atom->dvector[isite1][id]/nTotal;
+ nMoleculesOld1 = atom->dvector[isite1+nspecies][id];
+ nMolecules1 = atom->dvector[isite1][id];
+ fractionOld1 = nMoleculesOld1/nTotalOld;
+ fraction1 = nMolecules1/nTotal;
}
if (isOneFluid(isite2) == false){
- fractionOld2 = atom->dvector[isite2+nspecies][id]/nTotalOld;
- fraction2 = atom->dvector[isite2][id]/nTotal;
+ nMoleculesOld2 = atom->dvector[isite2+nspecies][id];
+ nMolecules2 = atom->dvector[isite2][id];
+ fractionOld2 = nMoleculesOld2/nTotalOld;
+ fraction2 = nMolecules2/nTotal;
}
if (isOneFluid(isite1) || isOneFluid(isite2)){
- fractionOld = 0.0;
- fraction = 0.0;
+ nMoleculesOFAold = 0.0;
+ nMoleculesOFA = 0.0;
+ fractionOFAold = 0.0;
+ fractionOFA = 0.0;
for (int ispecies = 0; ispecies < nspecies; ispecies++){
if (isite1 == ispecies || isite2 == ispecies) continue;
- fractionOld += atom->dvector[ispecies+nspecies][id] / nTotalOld;
- fraction += atom->dvector[ispecies][id] / nTotal;
+ nMoleculesOFAold += atom->dvector[ispecies+nspecies][id];
+ nMoleculesOFA += atom->dvector[ispecies][id];
+ fractionOFAold += atom->dvector[ispecies+nspecies][id] / nTotalOld;
+ fractionOFA += atom->dvector[ispecies][id] / nTotal;
}
if (isOneFluid(isite1)){
- fractionOld1 = fractionOld;
- fraction1 = fraction;
+ nMoleculesOld1 = 1.0-(nTotalOld-nMoleculesOFAold);
+ nMolecules1 = 1.0-(nTotal-nMoleculesOFA);
+ fractionOld1 = fractionOFAold;
+ fraction1 = fractionOFA;
}
if (isOneFluid(isite2)){
- fractionOld2 = fractionOld;
- fraction2 = fraction;
+ nMoleculesOld2 = 1.0-(nTotalOld-nMoleculesOFAold);
+ nMolecules2 = 1.0-(nTotal-nMoleculesOFA);
+ fractionOld2 = fractionOFAold;
+ fraction2 = fractionOFA;
}
}
+
+ if(fractionalWeighting){
+ mixWtSite1old = fractionOld1;
+ mixWtSite1 = fraction1;
+ mixWtSite2old = fractionOld2;
+ mixWtSite2 = fraction2;
+ } else {
+ mixWtSite1old = nMoleculesOld1;
+ mixWtSite1 = nMolecules1;
+ mixWtSite2old = nMoleculesOld2;
+ mixWtSite2 = nMolecules2;
+ }
}
/* ---------------------------------------------------------------------- */
diff --git a/src/USER-DPD/pair_multi_lucy_rx.h b/src/USER-DPD/pair_multi_lucy_rx.h
index 2913716c5a5ddce63670b869ff8064cc78ea78a3..87d2ed7ae8194f9fa0f03dfd85fdc0312461458d 100644
--- a/src/USER-DPD/pair_multi_lucy_rx.h
+++ b/src/USER-DPD/pair_multi_lucy_rx.h
@@ -78,7 +78,8 @@ class PairMultiLucyRX : public Pair {
int nspecies;
char *site1, *site2;
int isite1, isite2;
- void getParams(int, double &, double &, double &, double &);
+ void getMixingWeights(int, double &, double &, double &, double &);
+ bool fractionalWeighting;
};
diff --git a/src/USER-DPD/pair_table_rx.cpp b/src/USER-DPD/pair_table_rx.cpp
index 902d0e5bb4397144599637f8c03c977b7784ba6c..8498d752e31396b79101c362568beb06659d8b24 100644
--- a/src/USER-DPD/pair_table_rx.cpp
+++ b/src/USER-DPD/pair_table_rx.cpp
@@ -35,6 +35,12 @@ enum{NONE,RLINEAR,RSQ,BMP};
#define MAXLINE 1024
+#ifdef DBL_EPSILON
+ #define MY_EPSILON (10.0*DBL_EPSILON)
+#else
+ #define MY_EPSILON (10.0*2.220446049250313e-16)
+#endif
+
#define OneFluidValue (-1)
#define isOneFluid(_site_) ( (_site_) == OneFluidValue )
@@ -44,6 +50,7 @@ PairTableRX::PairTableRX(LAMMPS *lmp) : Pair(lmp)
{
ntables = 0;
tables = NULL;
+ fractionalWeighting = true;
}
/* ---------------------------------------------------------------------- */
@@ -82,21 +89,6 @@ void PairTableRX::compute(int eflag, int vflag)
if (eflag || vflag) ev_setup(eflag,vflag);
else evflag = vflag_fdotr = 0;
- double *fractionOld1, *fractionOld2;
- double *fraction1, *fraction2;
-
- {
- const int ntotal = atom->nlocal + atom->nghost;
-
- memory->create(fractionOld1, ntotal, "PairTableRx::compute::fractionOld1");
- memory->create(fractionOld2, ntotal, "PairTableRx::compute::fractionOld2");
- memory->create(fraction1, ntotal, "PairTableRx::compute::fraction1");
- memory->create(fraction2, ntotal, "PairTableRx::compute::fraction2");
-
- for (int i = 0; i < ntotal; ++i)
- getParams(i, fractionOld1[i], fractionOld2[i], fraction1[i], fraction2[i]);
- }
-
double **x = atom->x;
double **f = atom->f;
int *type = atom->type;
@@ -104,13 +96,29 @@ void PairTableRX::compute(int eflag, int vflag)
double *special_lj = force->special_lj;
int newton_pair = force->newton_pair;
- double fractionOld1_i, fractionOld1_j;
- double fractionOld2_i, fractionOld2_j;
- double fraction1_i, fraction1_j;
- double fraction2_i, fraction2_j;
+ double mixWtSite1old_i, mixWtSite1old_j;
+ double mixWtSite2old_i, mixWtSite2old_j;
+ double mixWtSite1_i, mixWtSite1_j;
+ double mixWtSite2_i, mixWtSite2_j;
double *uCG = atom->uCG;
double *uCGnew = atom->uCGnew;
+ double *mixWtSite1old = NULL;
+ double *mixWtSite2old = NULL;
+ double *mixWtSite1 = NULL;
+ double *mixWtSite2 = NULL;
+
+ {
+ const int ntotal = atom->nlocal + atom->nghost;
+ memory->create(mixWtSite1old, ntotal, "PairTableRx::compute::mixWtSite1old");
+ memory->create(mixWtSite2old, ntotal, "PairTableRx::compute::mixWtSite2old");
+ memory->create(mixWtSite1, ntotal, "PairTableRx::compute::mixWtSite1");
+ memory->create(mixWtSite2, ntotal, "PairTableRx::compute::mixWtSite2");
+
+ for (int i = 0; i < ntotal; ++i)
+ getMixingWeights(i, mixWtSite1old[i], mixWtSite2old[i], mixWtSite1[i], mixWtSite2[i]);
+ }
+
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
@@ -130,10 +138,10 @@ void PairTableRX::compute(int eflag, int vflag)
double uCGnew_i = 0.0;
double fx_i = 0.0, fy_i = 0.0, fz_i = 0.0;
- fractionOld1_i = fractionOld1[i];
- fractionOld2_i = fractionOld2[i];
- fraction1_i = fraction1[i];
- fraction2_i = fraction2[i];
+ mixWtSite1old_i = mixWtSite1old[i];
+ mixWtSite2old_i = mixWtSite2old[i];
+ mixWtSite1_i = mixWtSite1[i];
+ mixWtSite2_i = mixWtSite2[i];
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
@@ -147,10 +155,10 @@ void PairTableRX::compute(int eflag, int vflag)
jtype = type[j];
if (rsq < cutsq[itype][jtype]) {
- fractionOld1_j = fractionOld1[j];
- fractionOld2_j = fractionOld2[j];
- fraction1_j = fraction1[j];
- fraction2_j = fraction2[j];
+ mixWtSite1old_j = mixWtSite1old[j];
+ mixWtSite2old_j = mixWtSite2old[j];
+ mixWtSite1_j = mixWtSite1[j];
+ mixWtSite2_j = mixWtSite2[j];
tb = &tables[tabindex[itype][jtype]];
if (rsq < tb->innersq)
@@ -186,8 +194,8 @@ void PairTableRX::compute(int eflag, int vflag)
value = tb->f[itable] + fraction*tb->df[itable];
fpair = factor_lj * value;
}
- if (isite1 == isite2) fpair = sqrt(fractionOld1_i*fractionOld2_j)*fpair;
- else fpair = (sqrt(fractionOld1_i*fractionOld2_j) + sqrt(fractionOld2_i*fractionOld1_j))*fpair;
+ if (isite1 == isite2) fpair = sqrt(mixWtSite1old_i*mixWtSite2old_j)*fpair;
+ else fpair = (sqrt(mixWtSite1old_i*mixWtSite2old_j) + sqrt(mixWtSite2old_i*mixWtSite1old_j))*fpair;
fx_i += delx*fpair;
fy_i += dely*fpair;
@@ -208,11 +216,11 @@ void PairTableRX::compute(int eflag, int vflag)
((a*a*a-a)*tb->e2[itable] + (b*b*b-b)*tb->e2[itable+1]) *
tb->deltasq6;
if (isite1 == isite2){
- evdwlOld = sqrt(fractionOld1_i*fractionOld2_j)*evdwl;
- evdwl = sqrt(fraction1_i*fraction2_j)*evdwl;
+ evdwlOld = sqrt(mixWtSite1old_i*mixWtSite2old_j)*evdwl;
+ evdwl = sqrt(mixWtSite1_i*mixWtSite2_j)*evdwl;
} else {
- evdwlOld = (sqrt(fractionOld1_i*fractionOld2_j) + sqrt(fractionOld2_i*fractionOld1_j))*evdwl;
- evdwl = (sqrt(fraction1_i*fraction2_j) + sqrt(fraction2_i*fraction1_j))*evdwl;
+ evdwlOld = (sqrt(mixWtSite1old_i*mixWtSite2old_j) + sqrt(mixWtSite2old_i*mixWtSite1old_j))*evdwl;
+ evdwl = (sqrt(mixWtSite1_i*mixWtSite2_j) + sqrt(mixWtSite2_i*mixWtSite1_j))*evdwl;
}
evdwlOld *= factor_lj;
evdwl *= factor_lj;
@@ -238,10 +246,10 @@ void PairTableRX::compute(int eflag, int vflag)
}
if (vflag_fdotr) virial_fdotr_compute();
- memory->destroy(fractionOld1);
- memory->destroy(fractionOld2);
- memory->destroy(fraction1);
- memory->destroy(fraction2);
+ memory->destroy(mixWtSite1old);
+ memory->destroy(mixWtSite2old);
+ memory->destroy(mixWtSite1);
+ memory->destroy(mixWtSite2);
}
/* ----------------------------------------------------------------------
@@ -291,6 +299,8 @@ void PairTableRX::settings(int narg, char **arg)
else if (strcmp(arg[iarg],"msm") == 0) msmflag = 1;
else if (strcmp(arg[iarg],"dispersion") == 0) dispersionflag = 1;
else if (strcmp(arg[iarg],"tip4p") == 0) tip4pflag = 1;
+ else if (strcmp(arg[iarg],"fractional") == 0) fractionalWeighting = true;
+ else if (strcmp(arg[iarg],"molecular") == 0) fractionalWeighting = false;
else error->all(FLERR,"Illegal pair_style command");
iarg++;
}
@@ -1059,17 +1069,17 @@ double PairTableRX::single(int i, int j, int itype, int jtype, double rsq,
int tlm1 = tablength - 1;
Table *tb = &tables[tabindex[itype][jtype]];
- double fraction1_i, fraction1_j;
- double fraction2_i, fraction2_j;
- double fractionOld1_i, fractionOld1_j;
- double fractionOld2_i, fractionOld2_j;
+ double mixWtSite1_i, mixWtSite1_j;
+ double mixWtSite2_i, mixWtSite2_j;
+ double mixWtSite1old_i, mixWtSite1old_j;
+ double mixWtSite2old_i, mixWtSite2old_j;
fraction = 0.0;
a = 0.0;
b = 0.0;
- getParams(i,fractionOld1_i,fractionOld2_i,fraction1_i,fraction2_i);
- getParams(j,fractionOld1_j,fractionOld2_j,fraction1_j,fraction2_j);
+ getMixingWeights(i,mixWtSite1old_i,mixWtSite2old_i,mixWtSite1_i,mixWtSite2_i);
+ getMixingWeights(j,mixWtSite1old_j,mixWtSite2old_j,mixWtSite1_j,mixWtSite2_j);
if (rsq < tb->innersq) error->one(FLERR,"Pair distance < table inner cutoff");
@@ -1102,8 +1112,8 @@ double PairTableRX::single(int i, int j, int itype, int jtype, double rsq,
fforce = factor_lj * value;
}
- if (isite1 == isite2) fforce = sqrt(fraction1_i*fraction2_j)*fforce;
- else fforce = (sqrt(fraction1_i*fraction2_j) + sqrt(fraction2_i*fraction1_j))*fforce;
+ if (isite1 == isite2) fforce = sqrt(mixWtSite1_i*mixWtSite2_j)*fforce;
+ else fforce = (sqrt(mixWtSite1_i*mixWtSite2_j) + sqrt(mixWtSite2_i*mixWtSite1_j))*fforce;
if (tabstyle == LOOKUP)
phi = tb->e[itable];
@@ -1113,8 +1123,8 @@ double PairTableRX::single(int i, int j, int itype, int jtype, double rsq,
phi = a * tb->e[itable] + b * tb->e[itable+1] +
((a*a*a-a)*tb->e2[itable] + (b*b*b-b)*tb->e2[itable+1]) * tb->deltasq6;
- if (isite1 == isite2) phi = sqrt(fraction1_i*fraction2_j)*phi;
- else phi = (sqrt(fraction1_i*fraction2_j) + sqrt(fraction2_i*fraction1_j))*phi;
+ if (isite1 == isite2) phi = sqrt(mixWtSite1_i*mixWtSite2_j)*phi;
+ else phi = (sqrt(mixWtSite1_i*mixWtSite2_j) + sqrt(mixWtSite2_i*mixWtSite1_j))*phi;
return factor_lj*phi;
}
@@ -1141,46 +1151,74 @@ void *PairTableRX::extract(const char *str, int &dim)
/* ---------------------------------------------------------------------- */
-void PairTableRX::getParams(int id, double &fractionOld1, double &fractionOld2, double &fraction1, double &fraction2)
+void PairTableRX::getMixingWeights(int id, double &mixWtSite1old, double &mixWtSite2old, double &mixWtSite1, double &mixWtSite2)
{
- double nTotal = 0.0;
- double nTotalOld = 0.0;
+ double fractionOFAold, fractionOFA;
+ double fractionOld1, fraction1;
+ double fractionOld2, fraction2;
+ double nMoleculesOFAold, nMoleculesOFA;
+ double nMoleculesOld1, nMolecules1;
+ double nMoleculesOld2, nMolecules2;
+ double nTotal, nTotalOld;
+
+ nTotal = 0.0;
+ nTotalOld = 0.0;
for (int ispecies = 0; ispecies < nspecies; ++ispecies){
nTotal += atom->dvector[ispecies][id];
nTotalOld += atom->dvector[ispecies+nspecies][id];
}
- if(nTotal < 1e-8 || nTotalOld < 1e-8)
- error->all(FLERR,"The number of molecules in CG particle is less than 1e-8.");
+ if(nTotal < MY_EPSILON || nTotalOld < MY_EPSILON)
+ error->all(FLERR,"The number of molecules in CG particle is less than 10*DBL_EPSILON.");
if (isOneFluid(isite1) == false){
- fractionOld1 = atom->dvector[isite1+nspecies][id]/nTotalOld;
- fraction1 = atom->dvector[isite1][id]/nTotal;
+ nMoleculesOld1 = atom->dvector[isite1+nspecies][id];
+ nMolecules1 = atom->dvector[isite1][id];
+ fractionOld1 = nMoleculesOld1/nTotalOld;
+ fraction1 = nMolecules1/nTotal;
}
if (isOneFluid(isite2) == false){
- fractionOld2 = atom->dvector[isite2+nspecies][id]/nTotalOld;
- fraction2 = atom->dvector[isite2][id]/nTotal;
+ nMoleculesOld2 = atom->dvector[isite2+nspecies][id];
+ nMolecules2 = atom->dvector[isite2][id];
+ fractionOld2 = nMoleculesOld2/nTotalOld;
+ fraction2 = nMolecules2/nTotal;
}
if (isOneFluid(isite1) || isOneFluid(isite2)){
- double fractionOld = 0.0;
- double fraction = 0.0;
+ nMoleculesOFAold = 0.0;
+ nMoleculesOFA = 0.0;
+ fractionOFAold = 0.0;
+ fractionOFA = 0.0;
for (int ispecies = 0; ispecies < nspecies; ispecies++){
if (isite1 == ispecies || isite2 == ispecies) continue;
-
- fractionOld += atom->dvector[ispecies+nspecies][id]/nTotalOld;
- fraction += atom->dvector[ispecies][id]/nTotal;
+ nMoleculesOFAold += atom->dvector[ispecies+nspecies][id];
+ nMoleculesOFA += atom->dvector[ispecies][id];
+ fractionOFAold += atom->dvector[ispecies+nspecies][id]/nTotalOld;
+ fractionOFA += atom->dvector[ispecies][id]/nTotal;
}
-
if(isOneFluid(isite1)){
- fractionOld1 = fractionOld;
- fraction1 = fraction;
+ nMoleculesOld1 = 1.0-(nTotalOld-nMoleculesOFAold);
+ nMolecules1 = 1.0-(nTotal-nMoleculesOFA);
+ fractionOld1 = fractionOFAold;
+ fraction1 = fractionOFA;
}
-
if(isOneFluid(isite2)){
- fractionOld2 = fractionOld;
- fraction2 = fraction;
+ nMoleculesOld2 = 1.0-(nTotalOld-nMoleculesOFAold);
+ nMolecules2 = 1.0-(nTotal-nMoleculesOFA);
+ fractionOld2 = fractionOFAold;
+ fraction2 = fractionOFA;
}
}
-}
+ if(fractionalWeighting){
+ mixWtSite1old = fractionOld1;
+ mixWtSite1 = fraction1;
+ mixWtSite2old = fractionOld2;
+ mixWtSite2 = fraction2;
+ } else {
+ mixWtSite1old = nMoleculesOld1;
+ mixWtSite1 = nMolecules1;
+ mixWtSite2old = nMoleculesOld2;
+ mixWtSite2 = nMolecules2;
+ }
+}
diff --git a/src/USER-DPD/pair_table_rx.h b/src/USER-DPD/pair_table_rx.h
index f04ebced207fba9ad3cb8bd1e4b041aa13c367f8..c6afe6a8d59eea99f021add838112be13250683d 100644
--- a/src/USER-DPD/pair_table_rx.h
+++ b/src/USER-DPD/pair_table_rx.h
@@ -72,7 +72,8 @@ class PairTableRX : public Pair {
int nspecies;
char *site1, *site2;
int isite1, isite2;
- void getParams(int, double &, double &, double &, double &);
+ void getMixingWeights(int, double &, double &, double &, double &);
+ bool fractionalWeighting;
};
@@ -163,7 +164,7 @@ When using pair style table with a long-range KSpace solver, the
cutoffs for all atom type pairs must all be the same, since the
long-range solver starts at that cutoff.
-E: The number of molecules in CG particle is less than 1e-8
+E: The number of molecules in CG particle is less than 10*DBL_EPSILON
Self-explanatory. Check the species concentrations have been properly set
and check the reaction kinetic solver parameters in fix rx to more for