/DFT_pchem module | toolbox |============================== DFT_pchem =================================== | | Includes the HS and ES components. | | The main job of this module is to compute excess chemical potentials, | both in bulk and inhomogeneous situations. | | This version uses a newer formulation of the Rosenfeld expansion for the HS | aspect (the version used by Melchionna et al). (The polarity error present | in earlier versions of the vector terms has been corrected.) | |---------------------------------------------------------------------------- 200 dict dup begin | DFT_pchem |--------------------------- build pchem ------------------------------------ | pchempdict | -- /build_pchem { /pchempdict name |-- globals channel [ /SIGMA /Z /rCL /rCR /EU ] { Nspecies /d array def } forall [ /EPS ] { Ngrid /d array def } forall [ /Cs /EXs_hs /EXs_es /refCs /refEXs_es /E0s /EXs ] { [ Nspecies { Ngrid /d array } repeat ] def } forall [ /Ds ] { [ Nmobile { Ngrid /d array } repeat ] def } forall [ /bathDs /Ghydr ] { Nmobile /d array def } forall /selfFac Nspecies /d array def /refMSAradius Ngrid /d array def _channel |-- permanent locals pchempdict begin [ /CA /CB /CC /CD /CE /CF /AL /BT /Ct ] { Ngrid /d array def } forall [ /PNLs /MuNLs ] { [ 6 { Ngrid /d array } repeat ] def } forall /dMuNLs [ [ 0 3 Ngrid /d array ] | 0 [ 1 2 Ngrid /d array ] | 1 [ 1 3 Ngrid /d array ] | 2 [ 2 1 Ngrid /d array ] | 3 [ 2 2 Ngrid /d array ] | 4 [ 2 3 Ngrid /d array ] | 5 [ 2 5 Ngrid /d array ] | 6 [ 3 0 Ngrid /d array ] | 7 [ 3 1 Ngrid /d array ] | 8 [ 3 2 Ngrid /d array ] | 9 [ 3 3 Ngrid /d array ] | 10 [ 3 4 Ngrid /d array ] | 11 [ 3 5 Ngrid /d array ] | 12 [ 4 3 Ngrid /d array ] | 13 [ 4 5 Ngrid /d array ] | 14 [ 5 2 Ngrid /d array ] | 15 [ 5 3 Ngrid /d array ] | 16 [ 5 4 Ngrid /d array ] | 17 [ 5 5 Ngrid /d array ] | 18 ] def [ /dmun /dmun_w /dmun_ww ] { Ngrid /d array def } forall end } bind def /make_bulkdict { /n name 10 dict dup begin [ /Cs /EXs_hs /EXs_es ] { n /d array def } forall /PNLs 6 /d array def /MuNLs 6 /d array def end } bind def |============================ bulk MSA =========================== /msabuf 100 /d array def |------------------ compute bulk HS excess potentials, geometric measures, and | pressure from PY theory: | bulkdict | -- /comp_bulk_HS { /bulkdict name /nhs bulkdict /Cs get length def /nn 0 def /X0 sumps def /nn 1 def /X1 sumps def /nn 2 def /X2 sumps 4.0 Pi mul mul def /nn 3 def /X3 sumps 4.0 Pi mul 3.0 div mul def /xPY 1.0 X3 sub def /F0 xPY ln neg def /F1 X2 xPY div def /F2 X1 xPY div X2 xPY div 2 pwr 8.0 Pi mul div add def /F3 X1 X2 mul xPY div X0 add xPY div X2 xPY div 3 pwr 12.0 Pi mul div add def 0 1 nhs 1 sub { /i name /radius SIGMA i get 0.5 mul def F0 F1 radius mul add F2 4.0 mul Pi mul radius 2 pwr mul add F3 4.0 mul Pi mul 3.0 div radius 3 pwr mul add bulkdict /EXs_hs get i put } for X0 X1 X2 X3 0.0 0.0 5 -1 0 { /i name bulkdict /PNLs get i put } for F0 F1 F2 F3 0.0 0.0 5 -1 0 { /i name bulkdict /MuNLs get i put } for F3 bulkdict /Pressure put } bind def /sumps { 0.0 0 1 nhs 1 sub { /k name SIGMA k get 0.5 mul nn pwr bulkdict /Cs get k get mul add } for Avogadro mul } bind def |------------------ compute bulk phase ES excess potentials and screening length | bulkdict | -- /comp_bulk_ES { /bulkdict name /nes bulkdict /Cs get length def /contflag bulkdict /contflag get def msabuf [ /p /sig /t1 /t2 /t3 /t4 ] { exch nes /d parcel 3 -1 roll name } forall pop bulkdict /Cs get p copy Avogadro mul pop /sigma SIGMA 0 nes getinterval def /z Z 0 nes getinterval def /epses bulkdict /EPS get def find_Gamma comp_msa_par /lambda e0 dup mul kT div 4 div Pi div epses div epsilon0 div def sigma t1 copy Gamma mul 1.0 add pop z bulkdict /EXs_es get copy 2 pwr Gamma mul t1 div z t2 copy 2.0 mul sigma t3 copy 2 pwr t4 copy eta mul sub t1 div t3 eta 3.0 div mul add eta mul sigma mul add lambda mul neg pop Gamma bulkdict /Gamma put } bind def |--------------- find Gamma /comp_msa_par { /Delta 0.0 sigma t1 copy 3 pwr p mul add Pi -6.0 div mul 1.0 add def /Omega 0.0 t1 sig t2 copy Gamma mul 1.0 add div add Pi 2.0 div Delta div mul 1.0 add def /eta 0.0 p t1 copy sigma mul z mul t2 div add Pi 2.0 div Delta div mul Omega div def } bind def /find_Gamma { /prefac e0 e0 mul kT div epses div epsilon0 div def 1.0 0.0 sigma extrema pop /sigmin name sigma sig copy sigmin div pop contflag not { 0.0 z t1 copy 2 pwr p mul add prefac mul sqrt /Gamma name } if brent solveGamma begin { /brenttolerance 1e-8 def /guessfactor 0.1 def Gamma 0.5 mul sigmin mul 1e-2 100.0 1 guessbracbrent } stopped brentfct exch end _brent { (\nbrent did not find root\n) toconsole toconsole stop } { pop } ifelse sigmin div /Gamma name true bulkdict /contflag put } bind def /solveGamma 50 dict dup begin /brentfct ( Gamma in ES) def /eval { pchem /Gamma name /Delta 0.0 sigma t1 copy 3 pwr p mul add Pi -6.0 div mul 1.0 add def /Omega 0.0 t1 sig t2 copy Gamma mul 1.0 add div add Pi 2.0 div Delta div mul 1.0 add def /eta 0.0 p t1 copy sigma mul z mul t2 div add Pi 2.0 div Delta div mul Omega div def 0.0 sigma t1 copy 2 pwr eta mul neg z add t2 div 2 pwr p mul add prefac mul sigmin dup mul mul 4.0 Gamma dup mul mul sub _pchem } bind def end def |-------------------------- compute Donnan potential ----------------------- | These functions are needed by the root finder, brent, to solve the | electroneutrality equation at a boundary just left of index Ridx /solve_DonnanL 100 dict dup begin /brentmess (\nbrent: VdonnanL) def /eval { /Vdonnan name SF 0 1 Nmobile 1 sub { /kspecies name E0s kspecies get /E0 name EXs kspecies get /EX name DonnanC kspecies get E0 0 get EX 0 get add E0 poreproper dup length 2 div get EX poreproper dup length 2 div get add sub Vdonnan Z kspecies get mul sub exp mul Z kspecies get mul add } for } bind def end def /solve_DonnanR 100 dict dup begin /brentmess (\nbrent: VdonnanR) def /eval { /Vdonnan name SF 0 1 Nmobile 1 sub { /kspecies name E0s kspecies get /E0 name EXs kspecies get /EX name DonnanC kspecies get E0 last get EX last get add E0 poreproper dup length 2 div get EX poreproper dup length 2 div get add sub Vdonnan Z kspecies get mul sub exp mul Z kspecies get mul add } for } bind def end def |================================ DFT ================================= |----------------------------- build_DFT ------------------------------ | builds DFT variables/tables that depend on particle diameters | pchempdict | maxlanes (maximal span, in nodes, of particle interaction) /build_DFT { /pchempdict name make_hs_stencils pchempdict /HS_Stencils put make_es_stencils pchempdict /ES_Stencils put channel make_dEXs _channel } bind def |------------------------------- make the storage scheme for dEXs | NOTE: maximal number of lanes is determined only here /make_dEXs { /dEXs [ 0 1 Nspecies 1 sub { /kspecies1 name [ 0 1 Nspecies 1 sub { /kspecies2 name [ HS_Stencils kspecies1 get length 2 div HS_Stencils kspecies2 get length 2 div add 2 mul 1 add ES_Stencils kspecies1 Nspecies mul kspecies2 add get length 2 copy lt { exch } if pop { Ngrid /d array } repeat ] } for ] } for ] def 0 dEXs { { length 2 copy lt { exch } if pop } forall } forall } bind def |---------------------------- make stencils for HS convolutions | -- | HS_Stencils /make_hs_stencils { |---- stencils are designed separately for each species [ | start accumulation of stencils 0 1 Nspecies 1 sub { /kspecies name /radius SIGMA kspecies get SystemLength div 2 div def | scaled radius /jmax SIGMA kspecies get Xslice div 2 div ceil /l ctype def |-- determine the subranges of the grid where the left (right) half of each | particle species spans jmax,..,1 grid intervals size_stencils |-- create HS stencil storage of the species | elements are ordered from - to + in shifts /stencil [ jmaxm neg 1 -1 { /j name /kj j abs 1 sub def kml kj get /f name kmr kj get f sub 1 add /n name make_hs_element } for /j 0 def 0 /f name Ngrid /n name make_hs_element 1 1 jmaxp { /j name /kj j 1 sub def kpl kj get /f name kpr kj get f sub 1 add /n name make_hs_element } for ] def |-- negative j -1 -1 jmaxm neg { dup /j name abs 1 sub /kj name j jmaxm add /kel name |-- make truncated and full h's 0.0 WA copy pop /f kml kj get def /n kmr kj get f sub 1 add def X f j add 1 add n getinterval WA f n getinterval copy X f n getinterval sub radius add /ht name /hm H f j add 1 add n getinterval def j abs jmaxm lt { /kj kj 1 add def /ff kml kj get def /nn kmr kj get ff sub 1 add def H ff j add 1 add nn getinterval WA ff nn getinterval copy pop } if carve_some_hs |-- w3 /kw 3 def Pi neg Ca copy pop X f j add 1 add n getinterval Cb copy X f n getinterval sub Cc copy 2 pwr neg radius 2 pwr add Pi mul pop Cb -2.0 Pi mul mul pop left_hs_weight |-- w2 /kw 2 def 0 Ca copy Cb copy pop 2.0 Pi mul radius mul Cc copy pop left_hs_weight |-- wv2 /kw 5 def 0 Ca copy pop 2.0 Pi mul Cb copy pop X f j add 1 add n getinterval Cc copy X f n getinterval sub 2.0 Pi mul mul pop left_hs_weight } for |-- positive j 1 1 jmaxp { dup /j name 1 sub /kj name j jmaxp add /kel name |-- make truncated and full h's 0.0 WA copy pop /f kpl kj get def /n kpr kj get f sub 1 add def X f n getinterval WA f n getinterval copy radius add X f j add 1 sub n getinterval sub /ht name /hm H f j add n getinterval def j jmaxp lt { /kj kj 1 add def /ff kpl kj get def /nn kpr kj get ff sub 1 add def H ff j add nn getinterval WA ff nn getinterval copy pop } if carve_some_hs |-- w3 /kw 3 def Pi neg Ca copy pop X f j add 1 sub n getinterval Cb copy X f n getinterval sub Cc copy 2 pwr neg radius 2 pwr add Pi mul pop Cb -2.0 Pi mul mul pop right_hs_weight |-- w2 /kw 2 def 0 Ca copy Cb copy pop 2.0 Pi mul radius mul Cc copy pop right_hs_weight |-- wv2 /kw 5 def 0 Ca copy pop 2.0 Pi mul Cb copy pop X f j add 1 sub n getinterval Cc copy X f n getinterval sub 2.0 Pi mul mul pop right_hs_weight } for | positive j |-- compute derived weights and unscale stencil { 2 get /w name w 3 get SystemLength 2 pwr mul pop w 2 get SystemLength mul pop w 5 get SystemLength mul pop /rsradius radius SystemLength mul def w 2 get w 0 get copy 4.0 Pi mul rsradius 2 pwr mul div pop w 2 get w 1 get copy 4.0 Pi mul rsradius mul div pop w 5 get w 4 get copy 4.0 Pi mul rsradius mul div pop } forall |---- end of species loop stencil } for ] | close accumulation of stencils } bind def |------------------ carve some dynamic arrays /carve_some_hs { /al AL f n getinterval def /bt BT f n getinterval def /t WB f n getinterval def /Ca CA f n getinterval def /Cb CB f n getinterval def /Cc CC f n getinterval def } bind def |------------------ eqns 26-29: left hand weights /left_hs_weight { Ca al copy -0.25 mul ht mul Cb t copy 3 div add ht mul Cc t copy -0.5 mul add ht mul ht mul pop Ca bt copy 3 div ht mul Cb t copy -0.5 mul add ht mul Cc add ht mul pop stencil kel get 2 get kw get f n getinterval al hm div sub pop stencil kel 1 add get 2 get kw get f n getinterval al add bt add pop } bind def |------------------ eqns 12-14: right hand weights /right_hs_weight { Ca al copy 0.25 mul ht mul Cb t copy 3 div add ht mul Cc t copy 0.5 mul add ht mul ht mul pop Ca bt copy 3 div ht mul Cb t copy 0.5 mul add ht mul Cc add ht mul pop stencil kel get 2 get kw get f n getinterval al hm div add pop stencil kel 1 sub get 2 get kw get f n getinterval bt add al sub pop } bind def |------------------ make element of stencil /make_hs_element { [ |-- i_carver [ f n /getinterval mkact ] mkact bind |-- ij_carver [ f j add n /getinterval mkact ] mkact bind |-- weight arrays for PNL stencils (NOTE: these same weights can serve | for the EX stencils: weights 0-3 are used as are, and weights 4 and 5 | are negated; this utilizes the even/odd symmetries of the Rosenfeld | weight functions) [ 6 { 0.0 Ngrid /d array copy } repeat ] ] } bind def |------------------ compute truncated intervals /trunc_H { /ij_1 find mkpass dup 0 get /fel name 1 get /nel name j 0 lt { X fel 1 add nel getinterval WA ij_1 copy X fel j sub 1 add nel getinterval sub radius add } { X fel j sub 1 add nel getinterval WA ij_1 copy radius add X fel j sub 1 sub nel getinterval sub } ifelse pop } bind def |----------------------------- size stencils | /make_xx_element | -- /size_stencils { |-- for channel geometry [ /kml /kpl /kmr /kpr ] { jmax /l array def } forall 1 kml 0 put Ngrid 1 sub kmr 0 put 0 kpl 0 put Ngrid 2 sub kpr 0 put /i Ngrid 1 sub def /j -1 def { i j add 0 lt { /jmaxm j neg 1 sub def exit } if X i get radius sub X i j add get gt { /i i 1 sub def } { i kmr j abs put /j j 1 sub def j abs jmax ge { /jmaxm jmax def exit } if } ifelse } loop /i kmr jmaxm 1 sub get def /j jmaxm 1 sub neg def { X i get radius sub X i j add get lt { /i i 1 sub def } { i 1 add kml j abs put /j j 1 add def j -1 gt { exit } if } ifelse } loop /i 0 def /j 1 def { i j add Ngrid 2 sub ge { /jmaxp j 1 sub def exit } if X i get radius add X i j add get lt { /i i 1 add def } { i kpl j put /j j 1 add def j jmax ge { /jmaxp jmax def exit } if } ifelse } loop /i kpl jmaxp 1 sub get def /j jmaxp 1 sub def { X i get radius add X i j add get gt { /i i 1 add def } { i 1 sub kpr j put /j j 1 sub def j 1 lt { exit } if } ifelse } loop } bind def |----------------------------- make stencils for ES convolutions | NOTE: the weights EXCLUDE EPS | | -- | ES_Stencils /make_es_stencils { |---- stencils are designed separately for each species combination i, j [ | start accumulation of stencils 0 1 Nspecies 1 sub { /ispecies name 0 1 Nspecies 1 sub { /jspecies name /radius SIGMA ispecies get SIGMA jspecies get add 0.5 mul MSAradius dup add add SystemLength div def | scaled combined es radius /jmax radius SystemLength mul Xslice div ceil /l ctype def |-- determine the subranges of the grid where the left (right) half of each | particle species spans jmax,..,1 grid intervals size_stencils |-- create ES stencil storage of the species | elements are ordered from - to + in shifts /stencil [ jmaxm neg 1 -1 { /j name /kj j abs 1 sub def kml kj get /f name kmr kj get f sub 1 add /n name make_es_element } for /j 0 def 1 /f name Ngrid 2 sub /n name make_es_element 1 1 jmaxp { /j name /kj j 1 sub def kpl kj get /f name kpr kj get f sub 1 add /n name make_es_element } for ] def |-- compute the weights of the stencil |-- negative j -1 -1 jmaxm neg { dup /j name abs 1 sub /kj name j jmaxm add /kel name |-- make truncated and full h's 0.0 WA copy pop /f kml kj get def /n kmr kj get f sub 1 add def X f j add 1 add n getinterval WA f n getinterval copy X f n getinterval sub radius add /ht name /hm H f j add 1 add n getinterval def j abs jmaxm lt { /kj kj 1 add def /ff kml kj get def /nn kmr kj get ff sub 1 add def H ff j add 1 add nn getinterval WA ff nn getinterval copy pop } if carve_some_es |-- weight 1.0 3.0 div Ca copy pop radius Cb copy X f j add 1 add n getinterval add X f n getinterval sub Cc copy 2 pwr Cd copy Cb mul Ca mul pop left_es_weight } for |-- positive j 1 1 jmaxp { dup /j name 1 sub /kj name j jmaxp add /kel name |-- make truncated and full h's 0.0 WA copy pop /f kpl kj get def /n kpr kj get f sub 1 add def X f n getinterval WA f n getinterval copy radius add X f j add 1 sub n getinterval sub /ht name /hm H f j add n getinterval def j jmaxp lt { /kj kj 1 add def /ff kpl kj get def /nn kpr kj get ff sub 1 add def H ff j add nn getinterval WA ff nn getinterval copy pop } if carve_some_es |-- weight -1.0 3.0 div Ca copy pop radius Cb copy X f j add 1 sub n getinterval sub X f n getinterval add Cc copy 2 pwr neg Cd copy Cb mul Ca mul pop right_es_weight } for | positive j |-- compute derived weights and unscale /Cfactor Z ispecies get Z jspecies get mul e0 dup mul mul Avogadro mul 8.0 epsilon0 mul SIGMA ispecies get 0.5 mul MSAradius add mul SIGMA jspecies get 0.5 mul MSAradius add mul div SystemLength 4 pwr mul kT div def stencil { 2 get Cfactor mul pop } forall |---- end of species i,j loops stencil } for } for ] | close accumulation of stencils } bind def |------------------ carve some dynamic arrays /carve_some_es { /al AL f n getinterval def /bt BT f n getinterval def /t WB f n getinterval def /Ca CA f n getinterval def /Cb CB f n getinterval def /Cc CC f n getinterval def /Cd CD f n getinterval def } bind def |------------------ es eqns 65-66: left hand weights /left_es_weight { Ca al copy 0.2 mul ht mul Cb t copy 0.25 mul sub ht mul Cc t copy 3.0 div add ht mul Cd t copy 0.5 mul sub ht mul ht mul pop Ca bt copy -0.25 mul ht mul Cb t copy 3.0 div add ht mul Cc t copy -0.5 mul add ht mul Cd add ht mul pop stencil kel get 2 get f n getinterval al hm div sub pop stencil kel 1 add get 2 get f n getinterval al add bt add pop } bind def |------------------ es eqns 52-53: right hand weights /right_es_weight { Ca al copy 0.2 mul ht mul Cb t copy 0.25 mul add ht mul Cc t copy 3.0 div add ht mul Cd t copy 0.5 mul add ht mul ht mul pop Ca bt copy 0.25 mul ht mul Cb t copy 3.0 div add ht mul Cc t copy 0.5 mul add ht mul Cd add ht mul pop stencil kel get 2 get f n getinterval al hm div add pop stencil kel 1 sub get 2 get f n getinterval bt add al sub pop } bind def |------------------ make element of stencil /make_es_element { [ |-- i_carver [ f n /getinterval mkact ] mkact bind |-- ij_carver [ f j add n /getinterval mkact ] mkact bind |-- weight array 0.0 Ngrid /d array copy ] } bind def |------------------- convolve /convolve_es { /t name stencil { /element name element 0 get /i_iv name element 1 get /ij_iv name element 2 get /w name EX_es i_iv t ij_iv WA ij_iv copy w i_iv mul sub pop } forall } bind def |------------------------ compute excess potentials ------------------------ | needs Cs /comp_EXs { comp_PNLs comp_MuNLs comp_EXs_hs comp_EXs_es sum_EXs_components } bind def |----------------------------- sum HS and ES components of EXs /sum_EXs_components { 0 1 Nspecies 1 sub { /kspecies name EXs_hs kspecies get EXs kspecies get copy EXs_es kspecies get add pop } for } bind def |----------------------------- compute HS excess potentials /comp_EXs_hs { 0 1 Nspecies 1 sub { /kspecies name EXs_hs kspecies get /EX_hs name kspecies comp_DFT_hs_k Lbulk /EXs_hs get kspecies get EX_hs 0 put Rbulk /EXs_hs get kspecies get EX_hs last put } for } bind def |-- on selected species, in specified buffer /comp_DFT_hs_k { HS_Stencils exch get /stencil name 0.0 EX_hs copy pop 0 1 5 { /knl name stencil { /element name element 0 get /i_iv name element 1 get /ij_iv name element 2 get knl get /w name EX_hs i_iv MuNLs knl get ij_iv WA ij_iv copy w i_iv mul knl 4 lt { add } { sub } ifelse pop | SEE: make_hs_element } forall } for EX_hs SystemLength mul pop } bind def |-------------------------------- compute ES excess potentials /comp_EXs_es { 0 1 Nspecies 1 sub { /ispecies name /EX_es EXs_es ispecies get def 0.0 EX_es interior copy pop 0 1 Nspecies 1 sub { /jspecies name /ijspecies ispecies Nspecies mul jspecies add def /stencil ES_Stencils ijspecies get def refCs jspecies get WB copy neg Cs jspecies get ispecies isconfined ispecies jspecies eq and { Ct copy 1.0 selfFac ispecies get sub mul } if add convolve_es } for EX_es interior EPS interior div refEXs_es ispecies get interior add pop Lbulk /EXs_es get ispecies get EX_es 0 put Rbulk /EXs_es get ispecies get EX_es last put } for } bind def |------------------- compute derivatives of excess potentials with regard to Us | (expects that EXs have already been computed; needs also Cs) /comp_dEXs { comp_dMuNLs 0.0 dEXs { { { copy } forall } forall } forall pop comp_dEXs_hs comp_dEXs_es 0 1 Nspecies 1 sub { /kspecies1 name dEXs kspecies1 get /dEX1 name 0 1 Nspecies 1 sub { /kspecies2 name dEX1 kspecies2 get /dEX2 name dEX2 length /nlanes name 0 1 nlanes 1 sub { /klane name klane nlanes 2 div sub /kofs name kofs 0 le { /ij 0 def /n Ngrid kofs add def /i kofs neg def } if kofs 0 gt { /ij kofs def /n Ngrid kofs sub def /i 0 def } if dEX2 klane get i n getinterval Cs kspecies2 get ij n getinterval mul pop } for } for } for } bind def |-------------------- fill in derivatives of HS potentials /old_comp_dEXs_hs { dMuNLs { /dMuNL name dMuNL 0 get /kmu name dMuNL 1 get /kn name dMuNL 2 get dmun copy Avogadro SystemLength 2 pwr mul mul pop 0 1 Nspecies 1 sub { /kspecies1 name dEXs kspecies1 get /dEX1 name HS_Stencils kspecies1 get /stencil1 name stencil1 length /nlanes1 name 0 1 nlanes1 1 sub { /klane1 name stencil1 klane1 get /lane1 name lane1 0 get /i_iv1 name lane1 1 get /ij_iv1 name lane1 2 get kmu get /w1 name 0.0 dmun_w copy pop dmun ij_iv1 dmun_w ij_iv1 copy w1 i_iv1 mul pop 0 1 Nspecies 1 sub { /kspecies2 name dEX1 kspecies2 get /dEX2 name HS_Stencils kspecies2 get /stencil2 name stencil2 length /nlanes2 name 0 1 nlanes2 1 sub { /klane2 name stencil2 klane2 get /lane2 name lane2 2 get kn get /w2 name w2 dmun_ww copy dmun_w mul pop dEX2 klane1 nlanes1 2 div sub klane2 nlanes2 2 div sub add dEX2 length 2 div add get i_iv1 dmun_ww ij_iv1 add pop } for } for } for } for } forall } bind def |-- computation of the HS derivatives, optimized to get rid of extraneous | operations; refer to non-optimized code above for clarity /comp_dEXs_hs { dMuNLs { /dMuNL name dMuNL 0 get /kmu name dMuNL 1 get /kn name dMuNL 2 get dmun copy Avogadro SystemLength 2 pwr mul mul pop 0 1 Nspecies 1 sub { /kspecies1 name dEXs kspecies1 get /dEX1 name HS_Stencils kspecies1 get /stencil1 name stencil1 length /nlanes1 name /half_nlanes1 nlanes1 2 div def 0 1 nlanes1 1 sub { /klane1 name /knlane1 klane1 half_nlanes1 sub def stencil1 klane1 get /lane1 name lane1 0 get /i_iv1 name lane1 1 get /ij_iv1 name lane1 2 get kmu get i_iv1 dmun_w ij_iv1 copy dmun ij_iv1 mul pop 0 1 Nspecies 1 sub { /kspecies2 name HS_Stencils kspecies2 get /stencil2 name stencil2 length /nlanes2 name dEX1 kspecies2 get dup /dEX2 name length 2 div knlane1 add nlanes2 2 div sub /idx_dEX2 name 0 1 nlanes2 1 sub { /klane2 name stencil2 klane2 get 2 get kn get ij_iv1 dmun_ww ij_iv1 copy dmun_w ij_iv1 mul dEX2 idx_dEX2 klane2 add get i_iv1 exch add pop } for } for } for } for } forall } bind def |---------------------- fill in derivatives of ES potentials /comp_dEXs_es { 0 1 Nspecies 1 sub { /kspecies1 name dEXs kspecies1 get /dEX1 name 0 1 Nspecies 1 sub { /kspecies2 name dEX1 kspecies2 get /dEX2 name ES_Stencils kspecies1 Nspecies mul kspecies2 add get /stencil name /nlanes stencil length def 0 1 nlanes 1 sub { /klane name stencil klane get /lane name lane 0 get /i_iv name lane 2 get /w name dEX2 klane nlanes 2 div sub dEX2 length 2 div add get i_iv w i_iv dmun i_iv copy EPS i_iv div kspecies1 isconfined kspecies1 kspecies2 eq and { 1.0 selfFac kspecies1 get sub mul } if sub pop } for } for } for } bind def |---------------------------- DFT support ----------------------------- |------------------- compute excess potentials for non-locals /comp_MuNLs { 1.0 CA copy PNLs 3 get sub pop | 1-n3 PNLs 5 get CD copy PNLs 2 get div dup mul neg 1.0 add 3 pwr pop | RT |-- mu_0 CA MuNLs 0 get copy ln neg pop CA -1 pwr pop | 1/(1-n3) |-- mu_1 PNLs 2 get MuNLs 1 get copy CA mul |-- mu_2 CC copy 2 pwr 8.0 Pi mul -1 pwr mul | CC MuNLs 2 get copy CD mul PNLs 1 get CE copy CA mul add CC PNLs 2 get mul 3.0 -1 pwr mul pop | CC PNLs 5 get CB copy PNLs 2 get div dup mul | CB CE copy CC mul 6.0 mul PNLs 2 get div CB neg 1.0 add 2 pwr mul add pop | CB |-- mu_3 CC MuNLs 3 get copy 2.0 mul CD mul PNLs 0 get add PNLs 1 get CE copy PNLs 2 get mul CA mul add PNLs 4 get CE copy PNLs 5 get mul CA mul sub CA mul pop |-- mu_4 PNLs 5 get MuNLs 4 get copy CA mul neg pop |-- mu_5 CC MuNLs 5 get copy -6.0 mul PNLs 5 get mul CB mul PNLs 2 get CB copy -2 pwr mul PNLs 4 get CC copy CA mul sub pop |-- plug in endpoints of grid 0 1 5 { /knl name Lbulk /MuNLs get knl get MuNLs knl get 0 put Rbulk /MuNLs get knl get MuNLs knl get last put } for } bind def |------------------- compute derivatives of non-local excess potentials /comp_dMuNLs { PNLs 3 get CA copy neg 1.0 add -1 pwr | CA = 1/(1-n3) dMuNLs 0 get 2 get copy dMuNLs 1 get 2 get copy dMuNLs 3 get 2 get copy dMuNLs 7 get 2 get copy dMuNLs 14 get 2 get copy dMuNLs 17 get 2 get copy pop CA CB copy dup mul pop | CB = 1/(1-n3)^2 CB CC copy PNLs 2 get mul pop | CC = n2/(1-n3)^2 CC dMuNLs 2 get 2 get copy dMuNLs 8 get 2 get copy pop CB CC copy PNLs 5 get mul neg | CC = -n5/(1-n3)^2 dMuNLs 11 get 2 get copy dMuNLs 13 get 2 get copy neg pop CB CC copy PNLs 1 get mul | CC = n1/(1-n3)^2 dMuNLs 5 get 2 get copy dMuNLs 9 get 2 get copy pop PNLs 5 get CB copy dup mul | CB = a CC copy PNLs 2 get CD copy dup mul sub pop | CC = a-n2^2 PNLs 2 get CD copy 2 pwr | CD = n2^2 dMuNLs 4 get 2 get copy CB add CD mul 2.0 CE copy CB mul CB mul add CC mul neg 1.0 CE copy PNLs 3 get sub CD mul 2 pwr PNLs 2 get mul 4.0 Pi mul mul div pop dMuNLs 5 get 2 get CC CE copy CD div CA mul 2 pwr CA mul 4.0 Pi mul -1 pwr mul CD CF copy CB add mul add pop 3.0 dMuNLs 6 get 2 get copy CB mul CD add PNLs 5 get mul CC mul 1.0 CE copy PNLs 3 get sub CD mul 2 pwr 4.0 Pi mul mul div pop dMuNLs 9 get 2 get PNLs 2 get CD copy 2 pwr CB add CC CE copy CA mul PNLs 2 get div PNLs 2 get div dup mul CA mul mul 4.0 Pi mul -1 pwr mul add pop CC dMuNLs 10 get 2 get copy PNLs 2 get div 3 pwr CA mul 4.0 Pi mul -1 pwr mul neg PNLs 1 get CD copy PNLs 2 get mul PNLs 4 get CE copy PNLs 5 get mul sub 2.0 mul add CA mul PNLs 0 get add CA mul CA mul pop CC dMuNLs 12 get 2 get copy PNLs 2 get div 2 pwr PNLs 2 get div CA mul PNLs 5 get mul -2.0 Pi mul -1 pwr mul PNLs 4 get sub CA mul CA mul pop CA dMuNLs 15 get 2 get copy PNLs 2 get CD copy -2 pwr mul 2 pwr PNLs 5 get mul 4.0 Pi mul -1 pwr mul CC mul PNLs 2 get CD copy 2 pwr CB CE copy 3.0 mul add mul neg pop CC dMuNLs 16 get 2 get copy PNLs 2 get div 2 pwr PNLs 2 get div CA mul PNLs 5 get mul -2.0 Pi mul -1 pwr mul PNLs 4 get sub CA mul CA mul neg pop PNLs 2 get dMuNLs 18 get 2 get copy -1 pwr CA mul 2 pwr PNLs 2 get div -4.0 Pi mul -1 pwr mul CC mul PNLs 2 get CD copy 2 pwr neg CB 5.0 mul add mul neg pop } bind def |------------------- compute non-local densities /comp_PNLs { 0.0 PNLs { copy } forall pop 0 1 Nspecies 1 sub { /kkspecies name Cs kkspecies get /C name HS_Stencils kkspecies get /stencil name stencil { /element name element 0 get /i_iv name element 1 get /ij_iv name element 2 get /ws name 0 1 5 { /kw name PNLs kw get i_iv C ij_iv WA ij_iv copy ws kw get i_iv mul add pop } for } forall } for |-- unscale and establish bulk values at ends of grid 0 1 5 { /kn name PNLs kn get /PNL name PNL Avogadro SystemLength mul mul pop Lbulk /PNLs get kn get PNL 0 put Rbulk /PNLs get kn get PNL last put } for } bind def |---------------------------------- tests ------------------------------------ /test_dEXs { /kknode2 name /kkspecies2 name /kknode1 name /kkspecies1 name /oofset kknode2 kknode1 sub def dEXs kspecies1 get kspecies2 get dup length 2 div oofset add get kknode1 get _ pop /temp Cs 0 get 0 get 1e-3 add def /dh temp Cs 0 get 0 get sub def Cs kspecies2 get kknode2 get dh add Cs kspecies2 get kknode2 put comp_EXs /valp EXs kspecies1 get kknode1 get def Cs kspecies2 get kknode2 get dh 2.0 mul sub Cs kspecies2 get kknode2 put comp_EXs /valm EXs kspecies1 get kknode1 get def Cs kspecies2 get kknode2 get dh add Cs kspecies2 get kknode2 put comp_EXs valp valm sub dh 2.0 mul div Cs kspecies2 get kknode2 get mul _ pop } bind def end _module