怎样用probit计算半数有效浓度
Posted
tags:
篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了怎样用probit计算半数有效浓度相关的知识，希望对你有一定的参考价值。
参考技术A 药物安全性指标 EC50：引起50%个体有效的剂量扩大多少倍，才会导致50%的个体死亡（或中毒）。此值越大，用药相对越安全。 在医药卫生科研等工作中，经常需要测定某种因素（如：药物、毒物、细菌等）对实验动物（或昆虫等）的毒力或效力。杀虫剂毒力测定中常用的指标有半数致死量（LD50）、半数致死浓度（LC50）、半数有效量（ED50）等，现在广泛采用机率值分析法计算。本回答被提问者采纳计算机解决高中离子浓度计算
利用简单的化学及编程知识解决一类水溶液中各粒子浓度的计算
化学理论部分
假设氢离子浓度
假设现在有 \\(n\\) 元酸——\\(\\ceH_nAc\\)，其每一步电离的平衡常数为 \\(K_a_1,K_a_2,...,K_a_n\\)。（尽管 \\(n\\) 是个非常小的数，还是写成这种较为通用的形式）
那么如果我们知道了当前溶液中的氢离子浓度 \\(c(\\ceH)\\)，和酸根离子（或酸）的浓度 \\(c(\\ceH_yAc^(n-y)-)\\)。可以得出以下结论。
以三元酸 \\(\\ceH3PO4\\) 举例。写出以下电离方程式
于是得出
那么所有存在形式的浓度均可以用 \\(c(\\cePO_4^3-)\\) 和 \\(c(\\ceH^+)\\)表示。如果我们用想用 \\(c(\\ceH3PO4)\\) 来表示，则可以得到 \\(c(\\cePO_4^3-)=K_a_1K_a_2K_a_3\\times c(\\ceH^+)^-3\\times c(\\ceH3PO4)\\)。
不妨设往溶液中加入的磷酸的浓度为 \\(c_总\\)，就用 \\(c(\\cePO4^3-)\\) 作为未知数，利用物料守恒可以得到方程
显然我们如果知道了 \\(c(\\ceH^+)\\)，和 \\(c_总\\)，可以轻易地计算出 \\(c(\\cePO4^3-)\\)，进而计算出其他离子的浓度。
那么我们就可以将四种含 \\(\\ceP\\) 粒子的浓度随 \\(p\\ce H\\) 变化的图像画出来，这就是在题目中常见的图。
如果溶液中还有其他的酸或碱，其各个离子的浓度也可以根据这种方法计算出来。
考虑电荷守恒
放在实际的溶液中，其电荷一定是守恒的。我们上方只得到了各粒子浓度关于氢离子浓度的函数，正负电荷并不一定守恒。于是我们现在需要计算出一个合适的氢离子浓度，使得溶液中电荷守恒，那么这种情况就是满足的。由于各粒子浓度均非负数，那么我们应该可以证明氢离子浓度的解是唯一的（我没有证明，感觉可能是唯一的）。
程序设计部分
如何计算出氢离子浓度
我们根据电荷守恒，可以得到一个一元高次方程。如果我们尝试整理系数并利用求根公式求解，那么在程序上可能会比较难以实现。
不妨尝试试出氢离子浓度，计算出各阴阳离子的浓度之后，直接比较等式左右电荷绝对值的差，如果这个差的绝对值小于一个数，我们可以认为这个氢离子浓度是正确的。
但是如果按照浓度逐一尝试，每次为尝试的氢离子浓度增加一定的浓度，那么如果我们要保证精度，就需要尝试非常多次。
可以感性理解出正负电荷量差的绝对值随氢离子浓度的变化是一个连贯的曲线，我们的目标是找到这个曲线上的最小值点，这个点函数值为 \\(0\\)。这个函数上面可能有无数极值点（还是我感觉的，也有可能是单调的，这个时候可以直接使用二分法解决，显然但粒子种类非常多的时候已经很难用手证明是否在整个定义域上都是单调的了），而我们需要找到那个最小值点。这个时候可以考虑使用模拟退火算法（退火是不是本身就是化学名词）。
模拟退火就是比较常规的写法，记得比较上一个位置的答案而不是目前已经得到的最优的答案（很多模拟退火的竞赛题都是直接比较最优的答案，但是这里不行）。
一个精度问题
还是以磷酸举例，如果我们以 \\(c(\\cePO4^3-)\\) 为中间变量计算，如果 \\(c(\\ceH^+)\\) 的浓度非常大，那么得到的 \\(c(\\cePO4^3-)\\) 会非常小，由于默认的小数存储结构精度有限，一些有效数字可能会被直接省略掉，造成较大的精度丢失。那么我们可以在计算之前先找出在一定氢离子浓度下占据主导地位的粒子，然后用直接用它作为中间变量计算。（也许这点误差真的能够忽略不计吧，但是还是在这里写出来吧）
另一个大问题
磷酸 \\(K_a_1=6.9\\times 10^-3,K_a_2=6.2\\times 10^-8,K_a_3=4.8\\times 10^-13\\)，乘起来是一个非常小的数字，超过了 C 语言所能记录的浮点数精度。这意味着在计算过程中极易出现爆精度的问题，尽管可以采用先乘后除的办法，并且这些极小的数所代表的粒子浓度本身就可以忽略不计，但是还是会在一些玄学的地方带来误差，这样在比较两个浓度都非常小的粒子的时候就不能保证准确性了（我瞎猜说）。
解决的办法就是写自己写记录小数的数据结构（即高精度），但是这会浪费非常多的时间，并转化为一个彻头彻尾的代码问题。
实际上，程序中采用的课本上的平衡常数本身就可能存在着误差，所以在这里我们并不考虑解决这个问题。
一些目前版本的bug
- 对于在实际情况中浓度相同的粒子可能计算出的结果会有细微差异。
- 在我的电脑上运行程序的时候载入会有几秒钟的时间（这确实有点不太正常），可能是代码输入部分的实现问题。我不是专业的程序员，完全不会写 IO，请见谅。
一些省流
\\(0.1\\cemol\\cdot L^-1\\ceH3PO4\\) 溶液的计算结果：
pH: 1.6375
c(H_3PO_4): 0.076956557415177
c(H^+): 0.023043504584859
c(H_2PO_4^-): 0.023043380585157
c(HPO_4^2-): 0.000000061999666
c(OH^-): 0.000000000000434
c(PO_4^3-): 0.000000000000000
\\(0.1\\cemol\\cdot L^-1\\ceCH3COONa\\) 溶液的计算结果：
pH: 8.8785
c(Na^+): 0.100000000000000
c(CH_3COO^-): 0.099992441663165
c(OH^-): 0.000007559670777
c(CH_3COOH): 0.000007558336835
c(H^+): 0.000000001322809
\\(0.1\\cemol\\cdot L^-1\\ceNH3\\cdot H2O\\) 溶液的计算结果：
pH: 11.1247
c(NH_3-H_2O): 0.098667174249750
c(NH_4^+): 0.001332825750250
c(OH^-): 0.001332514123592
c(H^+): 0.000000000007505
代码
建立结构体记录每种弱酸或碱。结构体内有根据氢离子浓度计算出对应酸根离子们的浓度的函数。
#include <bits/stdc++.h>
#include <iomanip>
#include <iostream>
using namespace std;
typedef long double ld;
typedef pair<string,ld>ttfa;//这个 pair 记录答案中的离子名称及相应的浓度
#define mp make_pair
const ld K_w=1e-14;
mt19937 gen(chrono::system_clock::now().time_since_epoch().count());
uniform_real_distribution<ld>rd(0,1);
struct basecid
int n;ld c;//n 元酸或碱，以及浓度 c
vector<string>name;//名字，按照电荷数的绝对值从小到大排
vector<ld>K;//电离平衡常数
bool is_acid;//是酸(1)，还是碱(0)，如果是碱直接用将氢离子浓度转化成氢氧根离子浓度用同样的方法计算。
inline ld getk(int loc,ld c_H)//具体原理已经讲过，系数*该粒子浓度=总离子浓度，其中系数和氢离子及选择的粒子有关，这里求出这个系数
ld val=1,bas=1;
for(int i=loc-1;i>=0;--i)
bas=bas*c_H/K[i];
val+=bas;
bas=1;
for(int i=loc+1;i<=n;++i)
bas=bas*K[i-1]/c_H;
val+=bas;
return val;
inline ld e_it(ld c_H,ld c_all)//计算总电子数，c_all 如果为大于 0 则忽略目前结构体内的浓度，而使用这个浓度
if(c_all<0)c_all=c;
if(!is_acid)c_H=K_w/c_H;
int key=0;ld minv=getk(0,c_H);
for(int i=1;i<=n;++i)
ld now=getk(i,c_H);
if(now<minv)minv=now,key=i;
ld mid=c_all/minv,bas=mid,val=mid*key;
for(int i=key-1;i>=0;--i)
bas=bas*c_H/K[i];
bas=mid;
for(int i=key+1;i<=n;++i)
bas=bas*K[i-1]/c_H;
val+=i*bas;
return is_acid?-val:val;
inline vector<ttfa> lis_it(ld c_H,ld c_all)//单独又写了一个函数，将对应氢离子浓度下的各离子的浓度得出
if(c_all<0)c_all=c;
if(!is_acid)c_H=K_w/c_H;
int key=0;ld minv=getk(0,c_H);
for(int i=1;i<=n;++i)
ld now=getk(i,c_H);
if(now<minv)minv=now,key=i;
ld mid=c_all/minv,bas=mid;
vector<ttfa>lis(n+1);
lis[key]=name[key],mid;
for(int i=key-1;i>=0;--i)
bas=bas*c_H/K[i];
lis[i]=name[i],bas;
bas=mid;
for(int i=key+1;i<=n;++i)
bas=bas*K[i-1]/c_H;
lis[i]=name[i],bas;
return lis;
book[100];
struct normal//在溶液中浓度不变的粒子
string name;
int e;ld c;
inline ld e_it(ld c_all)
if(c_all<0)c_all=c;
return e*c_all;
;
vector<normal>other;
vector<basecid>mix;
inline ld e_forall(double c_H)//计算出该氢离子浓度下的电荷数
ld val=0;
for(auto x:mix)val+=x.e_it(c_H,-1);
for(auto x:other)val+=x.e_it(-1);
return val+c_H-K_w/c_H;
inline ld SA()//模拟退火，要记录 lastans。没怎么调参
ld loc=1e-7,ans=abs(e_forall(loc)),lasans=ans;
for(ld T=1;T>=1e-14;T*=0.995)
ld nloc=max(loc+(rd(gen)*2-1.0)*T,(ld)1e-16);
ld nval=abs(e_forall(nloc)),der=nval-lasans;
//printf("%Lf %Lf %Lf ",nloc,nval,der);
if(der<0)lasans=ans=nval,loc=nloc;
else if(exp(-der/T)>rd(gen))lasans=nval,loc=nloc/*,printf("accept")*/;
//puts("");
//printf("%.15Lf %.15Lf\\n",loc,ans);
return loc;
inline void initset()//这个写法可能有些丑陋，可能降低了效率
book[1].n=1,book[1].is_acid=1;
book[1].K.push_back(1.75e-5);
book[1].name.push_back("CH_3COOH");
book[1].name.push_back("CH_3COO^-");
book[2].n=1,book[2].is_acid=0;
book[2].K.push_back(1.8e-5);
book[2].name.push_back("NH_3-H_2O");
book[2].name.push_back("NH_4^+");
book[3].n=2,book[3].is_acid=1;
book[3].K.push_back(4.5e-7);
book[3].K.push_back(4.7e-11);
book[3].name.push_back("H_2CO_3");
book[3].name.push_back("HCO_3^-");
book[3].name.push_back("CO_3^2-");
book[4].n=3,book[4].is_acid=1;
book[4].K.push_back(6.9e-3);
book[4].K.push_back(6.2e-8);
book[4].K.push_back(4.8e-13);
book[4].name.push_back("H_3PO_4");
book[4].name.push_back("H_2PO_4^-");
book[4].name.push_back("HPO_4^2-");
book[4].name.push_back("PO_4^3-");
book[5].n=2,book[5].is_acid=1;
book[5].K.push_back(1.4e-2);
book[5].K.push_back(6.0e-8);
book[5].name.push_back("H_2SO_3");
book[5].name.push_back("HSO_3^-");
book[5].name.push_back("SO_3^2-");
book[6].n=1,book[6].is_acid=1;
book[6].K.push_back(4.0e-8);
book[6].name.push_back("HClO");
book[6].name.push_back("ClO-");
inline bool cmp(ttfa x,ttfa y)//这个写法比较丑陋，但是我并不会比较高级的写法
return x.second==y.second?x.first<y.first:x.second>y.second;
void forehead()
cout<<"English is used to avoid compatibility issues! (1.0 by BSE)"<<endl;
cout<<"----------------------------------------------------------------------------"<<endl;
cout<<"According to the design principle, this tool does not support direct input of compounds in the near future. However, the tool contain some normal particle, each one has a number."<<endl;
for(int i=1;i<=6;++i)
cout<<i<<": "<<book[i].name[0]<<endl;
cout<<"Note that they do not represent the concentration of the acid or base itself added to water, but rather the sum of the concentrations of all particles of this acid or base in the form in which they exist."<<endl;
cout<<"----------------------------------------------------------------------------"<<endl;
cout<<"If you want to add one of these particles, please enter the word auto first, then the number of it, at last the concentration."<<endl;
cout<<"If there are also some metal cations or acid ions of strong acids in the solution, please enter the word fixed first, then the name of it, followed by the charge number. At last the concentration."<<endl;
cout<<"If you want to add other weak acids or weak bases, please contact the writer! Because it\'s difficult to tell in English."<<endl;
cout<<"If you finish inputing, just enter the word end."<<endl;
cout<<"----------------------------------------------------------------------------"<<endl;
cout<<"For example, if you want to add 0.1mol/L NaCH_3COOH and 0.05mol/L NaCl, you just need to press: "<<endl;
cout<<"auto 1 0.1"<<endl<<"fixed Na^+ +1 0.15"<<endl<<"fixed Cl^- -1 0.05"<<endl<<"end"<<endl;
cout<<"----------------------------------------------------------------------------"<<endl;
void input()//没什么好说的，都是字面意思
while(1)
string opt;
cin>>opt;
if(opt=="auto")
int num;ld c;
cin>>num>>c;
basecid ac=book[num];
ac.c=c;
mix.push_back(ac);
else if(opt=="fixed")
normal fx;
cin>>fx.name>>fx.e>>fx.c;
other.push_back(fx);
else if(opt=="new")
basecid ac;
cin>>ac.n>>ac.is_acid;
for(int i=0;i<ac.n;++i)
ld K;cin>>K;
ac.K.push_back(K);
for(int i=0;i<=ac.n;++i)
string tmp;cin>>tmp;
ac.name.push_back(tmp);
cin>>ac.c;
mix.push_back(ac);
else if(opt=="end")break;
else cout<<"Non-conforming input!"<<endl;
return;
inline void pp()
ld c_H=SA();
vector<ttfa>each;
each.push_back(mp("H^+",c_H));each.push_back(mp("OH^-",K_w/c_H));
for(auto ac:mix)
vector<ttfa>tmp=ac.lis_it(c_H,-1);
for(auto ppair:tmp)
each.push_back(mp(ppair.first,ppair.second));
for(auto li:other)//这一步虽然可以在输入的时候完成，但是还是在这里写比较清晰一些
each.push_back(mp(li.name,li.c));
sort(each.begin(),each.end(),cmp);
cout<<"pH: "<<fixed<<setprecision(4)<<-log10(c_H)<<endl;
cout<<setiosflags(ios::fixed)<<setprecision(16);
for(auto ppair:each)//又是非常丑陋的写法
string tmp="c("+ppair.first+"): ";
cout<<setw(22)<<left<<tmp;
cout<<right<<ppair.second<<endl;
int main()
initset();
forehead();
input();cout<<"----------------------------------------------------------------------------"<<endl;
pp();
system("pause");
return 0;
后记
可能会改进一下算法，提升一下精度。最重要的是改良一下输入。
另外如果想要输入其他弱酸或弱碱，请按照下面的形式，看代码应该能知道怎么输入，下例子为 $ 0.1 \\cemol\\cdot L^-1\\ceH2C2O4$ 溶液。
new 2 1
5.6e-2 1.5e-4
H_2C_2O_4 HC_2O_4^- C_2O_4^2-
0.1
输出为
pH: 1.2836
c(H^+): 0.052049217002032
c(HC_2O_4^-): 0.051750936225208
c(H_2C_2O_4): 0.048099923386503
c(C_2O_4^2-): 0.000149140388288
c(OH^-): 0.000000000000192
先写这么多吧，还有好多想写的没写。如果错误和漏洞请多多反馈，谢谢观看。\\(\\rm by \\; BigSmall\\_En\\)
以上是关于怎样用probit计算半数有效浓度的主要内容，如果未能解决你的问题，请参考以下文章