科研科學(xué)計(jì)算高性能計(jì)算【免費(fèi)下載鏈接】lammpsPublic development project of the LAMMPS MD software package項(xiàng)目地址https://gitcode.com/gh_mirrors/la/lammps點(diǎn)擊查看免費(fèi)下載在含電解質(zhì)界面的介電系統(tǒng)中如疏水壁、膠體、脂質(zhì)雙分子層界面上會(huì)自發(fā)產(chǎn)生感應(yīng)電荷直接影響庫侖相互作用與離子力。LAMMPS 的fix polarize/bem/gmres、fix polarize/bem/icc、fix polarize/functional三個(gè) fix 命令隸屬 DIELECTRIC 包就是專門用于在兩個(gè)介電常數(shù)不同的介質(zhì)界面處計(jì)算并更新感應(yīng)電荷的工具。讀完本文你將掌握這三個(gè) fix 的語法與參數(shù)、atom_style dielectric數(shù)據(jù)文件的 15 列格式、fix_modify全部可調(diào)項(xiàng)、q_scaled輸出機(jī)制并能基于倉庫自帶的 confined / nopbc 示例直接搭建可運(yùn)行的界面極化模擬。一、命令語法與三種求解器三個(gè) fix 共用同一套語法ID與group-ID的含義與fix命令一致fix ID group-ID style nevery tolerance參數(shù)說明ID、group-IDfix 的標(biāo)識符與界面粒子所在原子組詳見 fix 命令文檔style取值為polarize/bem/gmres、polarize/bem/icc或polarize/functionalnevery該 fix 每隔多少個(gè)時(shí)間步被調(diào)用一次即ntimestep % nevery 0時(shí)才觸發(fā)求解源碼中pre_force有相應(yīng)判斷見 fix_polarize_bem_gmres.cpptolerance迭代求解器的相對收斂容差。文檔給出的三個(gè)最簡用法示例fix 2 interface polarize/bem/gmres 5 0.0001 fix 1 interface polarize/bem/icc 1 0.0001 fix 3 interface polarize/functional 1 0.0001三個(gè)求解器對應(yīng)三種獨(dú)立實(shí)現(xiàn)源碼分別位于 fix_polarize_bem_gmres.cpp、fix_polarize_bem_icc.cpp 和 fix_polarize_functional.cppDIELECTRIC 包 README 對其定位如下polarize/bem/gmres基于邊界元BEM方程、用 GMRES廣義最小殘差算法求解感應(yīng)電荷方法出自 Barros 等人的論文polarize/bem/icc同樣是邊界元方程但用逐次超松弛successive over-relaxation / SOR迭代求解方法出自 Tyagi 等人的論文polarize/functional直接對能量泛函做變分優(yōu)化方法出自 Jadhao 等人的論文。界面既可以是靜止的如接觸水的石墨烯層也可以是運(yùn)動(dòng)的如溶劑中的介電膠體——只要同時(shí)設(shè)置了合適的時(shí)間積分器例如fix rigid接口粒子會(huì)隨運(yùn)動(dòng)更新fix 每個(gè)nevery步重新求解感應(yīng)電荷。二、物理模型界面離散化與 dielectric 原子屬性這三個(gè) fix 計(jì)算的是兩不可穿透介質(zhì)界面處的感應(yīng)電荷。界面必須離散化為一系列頂點(diǎn)vertices每個(gè)頂點(diǎn)對應(yīng)一個(gè)邊界元且被當(dāng)作普通原子來存儲(chǔ)和通信。因此必須使用atom_style dielectric見 atom_style 文檔它額外定義了每個(gè)界面粒子的法向量、元素面積、局部介電失配等屬性——從 atom_vec_dielectric.h 可以看到該原子風(fēng)格新增的字段為mu法向量可通過mux/muy/muz像偶極矩分量一樣訪問、area、ed介電差、em介電均值、epsilon局部介電常數(shù)、curvature平均曲率和q_scaled。文檔用一個(gè)水-二氧化硅體系說明了各屬性的物理含義一側(cè)介質(zhì)介電常數(shù)為 78水另一側(cè)為 4二氧化硅界面離散為 2000 個(gè)邊界元每個(gè)界面試劑interface particle的法向量由二氧化硅指向水則法向介電差ed78 ? 4 74平均介電值em(78 4) / 2 41每個(gè)邊界元還帶有自己的面積和局部平均曲率fix 用它計(jì)算局部電場中的修正項(xiàng)若要模擬帶電界面界面試劑應(yīng)帶有非零電荷由其面積 × 表面電荷密度決定且其局部介電常數(shù)值取平均值em。對非界面粒子普通原子、帶電粒子而言法向量、面積、介電差字段無意義且不被使用它們的epsilon字段會(huì)被內(nèi)部使用計(jì)算庫侖相互作用時(shí)把真實(shí)電荷除以該值進(jìn)行縮放。例如模擬介電常數(shù)為 40 的隱式溶劑中攜帶 2 電荷的陽離子只需把q設(shè)為 2、epsilon設(shè)為 40無需手動(dòng)縮放電荷任何帶 dielectric 后綴的pair_style都會(huì)給出正確受力。文檔同時(shí)指出一個(gè)隱含假設(shè)由于epsilon不隨時(shí)間改變粒子在模擬中不能穿越界面。atom_style dielectric的 data 文件 Atoms 段包含15 列id mol type q x y z normx normy normz area/patch ed em epsilon curvature倉庫示例 data.confined 中一個(gè)界面行如下1 0 1 0 0 0 9.99798 0 0 1 0.866 8 6 6 0而離子行data.confined 末尾為4001 0 2 1 15 20 15 0 0 1 1.0 8 6 10 0 4002 0 3 -1 25 20 25 0 0 1 1.0 8 6 10 0按 示例 README 的說明normx/normy/normz是該邊界元界面處法單位向量的三個(gè)分量對真實(shí)離子無意義可填 0,0,1ed epsilon_out ? epsilon_in法向從epsilon_in指向epsilon_outem (epsilon_out epsilon_in)/2界面粒子的epsilon取em離子的epsilon取所在介質(zhì)的介電常數(shù)area/patch是元素面積對離子可填 1.0curvature為平均曲率球形界面取 1/半徑平面取 0。三、感應(yīng)電荷方程與三種算法3.1 邊界元方程對polarize/bem/gmres與polarize/bem/icc界面組內(nèi)頂點(diǎn)上的感應(yīng)電荷密度的表達(dá)式為$$\sigma_b(\mathbf{s}) \frac{1 - \bar{\epsilon}}{\bar{\epsilon}}, \sigma_f(\mathbf{s}) - \epsilon_0 \frac{\Delta \epsilon}{\bar{\epsilon}}, \mathbf{E}(\mathbf{s})\cdot\mathbf{n}(\mathbf{s})$$其中$\sigma_b$界面試劑 $\mathbf{s}$ 處的感應(yīng)電荷密度$\bar{\epsilon} (\epsilon_1 \epsilon_2)/2$界面試劑處平均介電常數(shù)即em$\Delta\epsilon \epsilon_1 - \epsilon_2$介電常數(shù)差即ed$\sigma_f$界面試劑處的自由電荷密度q/area$\mathbf{E}(\mathbf{s})$頂點(diǎn)處的電場$\mathbf{n}(\mathbf{s})$從 $\epsilon_2$ 介質(zhì)指向 $\epsilon_1$ 介質(zhì)的單位法向量。由于右端電場本身又依賴于感應(yīng)電荷這是一個(gè)隱式方程需要迭代求解。在源碼中右端向量 b 的構(gòu)造見 compute_induced_charges()其中每個(gè)頂點(diǎn)計(jì)算double ndotE epsilon0e2q * (Ex * norm[i][0] Ey * norm[i][1] Ez * norm[i][2]) / epsilon[i]; double sigma_f q[i] / area[i]; buffer[idx] (1 - em[i]) * sigma_f - ed[i] * ndotE / MY_4PI;epsilon0e2q是單位換算因子把 $\epsilon_0 \cdot E$ 換算為電荷/面積2單位源碼 setup() 中詳細(xì)注釋了real、metal、si、nano各單位制下該系數(shù)的推導(dǎo)并對不支持的單位制直接報(bào)錯(cuò)else if (strcmp(update-unit_style, lj) ! 0) error-all(FLERR, Only unit styles lj, real, metal, si and nano are supported);兩個(gè) BEM 求解器的迭代都在這兩類終止條件之一滿足時(shí)結(jié)束連續(xù)兩次迭代間感應(yīng)電荷的最大相對變化低于設(shè)定容差或迭代次數(shù)達(dá)到itr_max見下文fix_modify。polarize/bem/gmres采用 GMRES 算法Barros et al.在 gmres_solve() 中實(shí)現(xiàn)了 Arnoldi 正交化 Givens 旋轉(zhuǎn)的標(biāo)準(zhǔn)流程。值得注意的是實(shí)現(xiàn)上并不顯式存儲(chǔ)矩陣 A而是把作用于向量封裝為 apply_operator()把試探電荷寫入q_scaled調(diào)用 pair/kspace 風(fēng)格計(jì)算電場后讀回殘差則由 update_residual() 按 Barros 文獻(xiàn)式 (60) 直接更新。收斂判據(jù)在 gmres_solve() 尾部rho tol_rel * normbBarros 文獻(xiàn)建議的條件或rho rho_tol rho tol_abspolarize/bem/icc采用逐次超松弛算法Tyagi et al.。在 compute_induced_charges() 中每次迭代先用 pair/kspace 風(fēng)格重新計(jì)算電場然后做 SOR 更新q_bound (1 - omega) * q_bound omega * ((1.0 / em[i] - 1) * q_free - (ed[i] / em[i]) * ndotE * area[i]);收斂判據(jù)為各頂點(diǎn)感應(yīng)電荷相對變化的全局最大值MPI_Allreduce(..., MPI_MAX, ...)小于容差polarize/functional采用能量泛函變分法Jadhao et al.。從 源碼頭部注釋 看它求解的是線性系統(tǒng)(Rww Rww^T) w q Rwq并自 2019 年 10 月起切換為共軛梯度conjugate gradient求解器。該 fix 會(huì)在init()中向 neighbor 請求完整鄰居表NeighConst::REQ_FULLinit()其構(gòu)造函數(shù)的容差是可選參數(shù)缺省為 1.0e-6構(gòu)造函數(shù)。更完整的實(shí)現(xiàn)細(xì)節(jié)與推薦用法見 Nguyen 等人的論文(NguyenTD)Nguyen, Li, Bagchi, Solis, Olvera de la Cruz, Comput Phys Commun 241, 80-19 (2019)DIELECTRIC 包 README 也要求使用該包時(shí)引用該文獻(xiàn)。3.2 依賴的 pair / kspace 風(fēng)格文檔要求同時(shí)使用帶dielectric后綴的pair_style與kspace_style。源碼中 fix 通過dynamic_cast逐一識別兼容風(fēng)格update_efield() 列出的 pair 風(fēng)格有l(wèi)j/cut/coul/long/dielectric、lj/cut/coul/msm/dielectric、lj/cut/coul/cut/dielectric、lj/cut/coul/debye/dielectric、coul/long/dielectric、coul/cut/dielectric不匹配則報(bào) Pair style not compatiblekspace 風(fēng)格支持pppm/dielectric與msm/dielectric。電場計(jì)算按 pair 與 kspace 兩部分分別累積到efield_pair、efield_kspace再求和。主要庫侖風(fēng)格文檔見 pair lj/cut/coul/long/dielectric 與 kspace_style pppm/dielectric。四、q_scaled 輸出與 restart 行為三個(gè) fix 求出的感應(yīng)電荷存儲(chǔ)在q_scaled字段中。q字段始終是真實(shí)電荷、模擬中不改變對界面粒子q_scaled 真實(shí)電荷/局部介電常數(shù)epsilon 感應(yīng)電荷對非界面粒子q_scaled 真實(shí)電荷/epsilon。GMRES 實(shí)現(xiàn)中這一點(diǎn)體現(xiàn)得很直接compute_induced_charges() 尾部if (induced_charge_idx[i] 0) { int idx induced_charge_idx[i]; q_scaled[i] induced_charges[idx] * area[i] q[i]; } else { q_scaled[i] q_backup[i]; }注意感應(yīng)電荷induced_charges是電荷密度乘上area[i]才是該邊界元的電荷量。文檔給出的讀取示例compute qs all property/atom q_scaled dump 1 all custom 1000 all.txt id type q x y z c_qs關(guān)于 restart該 fix 的任何信息都不寫入二進(jìn)制 restart 文件感應(yīng)電荷會(huì)在setup()時(shí)重新求解初始化因此從 restart 恢復(fù)不影響正確性。polarize/bem/gmres與polarize/bem/icc還會(huì)計(jì)算一個(gè) 2 元素的全局向量global_freq 1size_vector 2構(gòu)造函數(shù)可通過各類輸出命令訪問第一個(gè)元素是求解器終止時(shí)的迭代次數(shù)上限為itr_max第二個(gè)元素是 RMS 誤差。compute_vector()的實(shí)現(xiàn)見 fix_polarize_bem_gmres.cpp#L969-L977可用variable x1 c_3[1]/c_3[2]3為 fix ID配合thermo_style custom或fix ave/time監(jiān)控收斂質(zhì)量。polarize/functional則不提供向量輸出。五、fix_modify 參數(shù)詳解通過 fix_modify 命令 可修改以下設(shè)置三個(gè) fix 的modify_param實(shí)現(xiàn)分別見 GMRES 與 ICCfix_modify ID itr_max arg fix_modify ID dielectrics ediff emean epsilon area charge fix_modify ID kspace arg fix_modify ID rand max seed fix_modify ID mr arg fix_modify ID omega argitr_max argarg 每步求解的最大迭代次數(shù)控制收斂迭代的上界dielectrics ediff emean epsilon area charge覆蓋group-ID中所有原子的介電屬性。ediff 介電差emean 介電均值epsilon 局部介電值area 元素面積charge 界面真實(shí)電荷任一參數(shù)傳NULL則保留 data 文件中的現(xiàn)有值。注意該設(shè)置作用于整個(gè)組——若兩個(gè)不同界面需要不同屬性需要為每個(gè)界面組分別建立一對fixfix_modify。源碼中NULL判斷與正數(shù)校驗(yàn)邏輯見 set_dielectric_params()kspace argyes或no開關(guān)長程相互作用參與感應(yīng)電荷計(jì)算默認(rèn) yes即若存在 kspace 風(fēng)格則自動(dòng)啟用若指定 yes 但系統(tǒng)沒有 kspace 風(fēng)格會(huì)打印警告 No Kspace pair style available 并禁用見 setup()rand max seed以均勻分布生成范圍 ?max/2 到max/2 的隨機(jī)初始感應(yīng)電荷隨機(jī)數(shù)種子為seed用于測試收斂性。源碼中實(shí)際會(huì)進(jìn)一步把隨機(jī)電荷全局減均值使總電荷為零init()mr arg僅polarize/bem/gmres求解表面電荷時(shí)使用的 q-向量Krylov 空間最大數(shù)量omega arg僅polarize/bem/iccTyagi 文獻(xiàn)定義的松弛參數(shù)一般應(yīng)設(shè)在 0 到 2 之間。默認(rèn)值匯總參數(shù)默認(rèn)值適用itr_max50全部kspaceyes全部omega0.7僅 ICCmr組內(nèi)原子數(shù) ? 1僅 GMRES隨機(jī)電荷初始化不啟用全部源碼印證itr_max 50在 GMRES 構(gòu)造函數(shù)、ICC 構(gòu)造函數(shù) 中初始化mr默認(rèn)在init()中被鉗制為mat_dim - 1fix_polarize_bem_gmres.cpp#L202-L203omega 0.7見 fix_polarize_bem_icc.cpp#L73。另外文檔提醒局部介電常數(shù)epsilon也可以獨(dú)立用 set 命令 設(shè)置示例腳本中正是這么做的。六、完整示例examples/PACKAGES/dielectric文檔指引的示例位于examples/PACKAGES/dielectric目錄包含兩套場景見 示例 READMEdata.confined兩個(gè)等量異號點(diǎn)電荷被限制在兩個(gè)平行界面之間epsilon12 | epsilon210 | epsilon12data.sphere兩個(gè)點(diǎn)電荷位于球形界面試劑之外epsilon_in1epsilon_out10。6.1 confined 示例周期性邊界 PPPM 長程入口腳本 in.confined.gmres 只有兩行通過變量選擇求解器后 include 公共腳本variable method index gmres # gmres BEM/GMRES # icc BEM/ICC* # dof Direct optimization of the functional include confined.incin.confined.icc 與 in.confined.dof 只需把method改為icc或dof。核心邏輯在 confined.inc完整流程如下# 界面法向量方向須與 ed 一致從 epsilon1 區(qū)指向 epsilon2 區(qū) variable epsilon1 index 20 variable epsilon2 index 10 variable data index data.confined newton off units lj atom_style dielectric atom_modify map array dimension 3 boundary p p f # 根據(jù)介電常數(shù)在腳本內(nèi)計(jì)算界面試劑屬性 variable ed equal v_epsilon2 - v_epsilon1 variable em equal (v_epsilon2 v_epsilon1)/2 variable epsilon equal 1.0 # 當(dāng)前未使用的 patch epsilon variable area equal 0.866 # patch 面積與 data 文件一致 read_data ${data} group interface type 1 group ions type 2 3 group cations type 2 group anions type 3 # 設(shè)置離子所在介質(zhì)的介電常數(shù)epsilon 用 set 命令獨(dú)立設(shè)置 set group cations epsilon ${epsilon2} set group anions epsilon ${epsilon2} pair_style lj/cut/coul/long/dielectric 1.122 10.0 pair_coeff * * 1.0 1.0 pair_coeff 1 1 0.0 1.0 kspace_style pppm/dielectric 0.0001 kspace_modify slab 3.0 neigh_modify every 1 delay 0 check yes one 5000 fix 1 ions nve # 用 fix / fix_modify 設(shè)置界面組屬性并按 method 選擇求解器 if ${method} gmres then fix 3 interface polarize/bem/gmres 1 1.0e-4 fix_modify 3 itr_max 50 dielectrics ${ed} ${em} ${epsilon} ${area} NULL elif ${method} icc fix 3 interface polarize/bem/icc 1 1.0e-4 fix_modify 3 itr_max 50 dielectrics ${ed} ${em} ${epsilon} ${area} NULL elif ${method} dof fix 3 interface polarize/functional 1 0.0001 fix_modify 3 dielectrics ${ed} ${em} ${epsilon} ${area} NULL else print Unsupported polarization solver thermo 1000 thermo_style custom step evdwl ecoul elong epair #f_3 thermo_modify flush yes run 0該腳本展示了文檔中提到的全部關(guān)鍵要素newton offfix 內(nèi)部假設(shè) pair 使用完整鄰居表并自行做通信、atom_modify map array配合 fix 內(nèi)按 atom tag 建立tag2mat映射、set命令設(shè)置離子epsilon、fix_modify dielectrics用變量覆蓋界面屬性。confined.inc 中還保留了被注釋掉的調(diào)試輸出可用compute efield/atom文檔導(dǎo)出離子與界面粒子上的電場分量來驗(yàn)證極化分布#compute ef all efield/atom #dump 1 all custom 100 all.dump id mol type q x y z #fx fy fz c_ef[1] c_ef[2] c_ef[3]6.2 nopbc 示例非周期邊界 大截?cái)鄋opbc.inc 演示無長程求解器的路線boundary f f fpair_style lj/cut/coul/cut/dielectric 1.122 20.0使用足夠大的截?cái)鄟砟M長程neigh_modify one 5000保證大截?cái)嘞锣従颖硗暾ＨN求解器的選擇方式與 confined 相同。README 給出的運(yùn)行方式為需已編譯 DIELECTRIC 包mpirun -np 8 ../../../src/lmp_mpi -in in.confined倉庫中同時(shí)保留著三種求解器分別的運(yùn)行日志如log.*.confined.gmres.g.4可直接對比同一場景下不同算法的收斂迭代數(shù)與能量。6.3 單元測試除示例外unittest/force-styles/tests/in.dielectric 配合 data.dielectric 被各 dielectric pair/kspace 風(fēng)格的 YAML 測試如 kspace-pppm_dielectric.yaml、mol-pair-lj_cut_coul_long_dielectric.yaml復(fù)用是驗(yàn)證atom_style dielectric與 dielectric 風(fēng)格受力一致性的回歸基準(zhǔn)修改介電相關(guān)代碼時(shí)可參考。七、構(gòu)建前提與使用限制7.1 包依賴這些 fix 屬于DIELECTRIC 包只有 LAMMPS 用該包編譯后才可用且必須同時(shí)安裝 KSPACE 包構(gòu)建細(xì)節(jié)見 Build package 文檔。依賴關(guān)系在 src/DIELECTRIC/Install.sh 中被硬校驗(yàn)# compute efield/atom is only usable when all styles # are installed, which in turn requires KSPACE if (test $1 1) then if (test ! -e ../pppm.cpp) then echo Must install KSPACE package with DIELECTRIC package exit 1 fi if (test ! -e ../pair_lj_cut_coul_debye.cpp) then echo Must install EXTRA-PAIR package with DIELECTRIC package exit 1 fi fi即安裝 DIELECTRIC 前系統(tǒng)里必須已有 KSPACEpppm.cpp與 EXTRA-PAIRpair_lj_cut_coul_debye.cpp文件。7.2 單位制與功能限制polarize/bem/gmres與polarize/bem/icc目前只支持lj、real、metal、si、nano五種 units源碼對其它單位制直接報(bào)錯(cuò)見上文 setup() 代碼polarize/functional尚不支持帶電界面界面粒子q非零的場景模型假設(shè)粒子不穿越界面因?yàn)閑psilon屬性在模擬中不改變兩個(gè) BEM fix 的fix_modify屬性覆蓋是組級別的同一 fix 管理的界面共享同一組ediff/emean/epsilon/area/charge。7.3 相關(guān)命令與參考資料相關(guān)命令pair_coeff、fix polarize、read_data、pair_style lj/cut/coul/long/dielectric、kspace_style pppm/dielectric、compute efield/atom輸出機(jī)制見 Howto output 與 restart。方法學(xué)參考文獻(xiàn)與文檔對應(yīng)(Barros)Barros, Sinkovits, Luijten, J. Chem. Phys, 140, 064903 (2014)——BEM 方程與 GMRES 求解(Tyagi)Tyagi, Suzen, Sega, Barbosa, Kantorovich, Holm, J Chem Phys, 132, 154112 (2010)——逐次超松弛ICC算法及omega參數(shù)(Jadhao)Jadhao, Solis, Olvera de la Cruz, J Chem Phys, 138, 054119 (2013)——能量泛函變分法(NguyenTD)Nguyen, Li, Bagchi, Solis, Olvera de la Cruz, Comput Phys Commun 241, 80-91 (2019)——本實(shí)現(xiàn)的完整描述與推薦用法使用 DIELECTRIC 包發(fā)表論文時(shí)應(yīng)引用。贊分享科研科學(xué)計(jì)算高性能計(jì)算【免費(fèi)下載鏈接】lammpsPublic development project of the LAMMPS MD software package項(xiàng)目地址https://gitcode.com/gh_mirrors/la/lammps點(diǎn)擊查看免費(fèi)下載相關(guān)推薦LAMMPS 中 compute efield/atom 命令詳解使用 DIELECTRIC 包計(jì)算原子級電場LAMMPS 中 compute efield/atom 命令詳解使用 DIELECTRIC 包計(jì)算原子級電場 本文以 LAMMPS 官方文檔 doc/src科研科學(xué)計(jì)算高性能計(jì)算LAMMPS compute dipole 命令詳解:體系偶極矩計(jì)算與 TIP4P 電荷位點(diǎn)修正LAMMPS compute dipole 命令詳解:體系偶極矩計(jì)算與 TIP4P 電荷位點(diǎn)修正 本文圍繞 LAMMPS 的 compute dipole 與科研科學(xué)計(jì)算高性能計(jì)算LAMMPS 介電常數(shù)設(shè)置指南dielectric 命令的語法、物理含義與源碼實(shí)現(xiàn)LAMMPS 介電常數(shù)設(shè)置指南dielectric 命令的語法、物理含義與源碼實(shí)現(xiàn) dielectric 是 LAMMPS 中用于全局調(diào)整庫侖Coulomb科研科學(xué)計(jì)算高性能計(jì)算上一篇GolangBypassAV 項(xiàng)目使用教程下一篇SignTools 開源項(xiàng)目教程創(chuàng)作聲明：本文部分內(nèi)容由AI輔助生成（AIGC），僅供參考