做了个新的反应活性位点的预测程序，欢迎大家试用

yjcmwgk

这半年多，做了个新的反应活性预测算法，带程序的。
我自认为它能准确预测亲电/亲核反应活性的具体位点和强度。

程序的exe文件在这儿下载
https://www.rsc.org/suppdata/d5/re/d5re00506j/d5re00506j1.zip

备注：有人说下载链接有时候打不开，那你直接去https://doi.org/10.1039/D5RE00506J，点“Supplementary files”里的zip文件，就是软件本身。
我喜欢把程序作为Supplementary files，和论文一起发表。
明面上的理由：为了让大家复现计算结果更方便！（这是真的）
暗戳戳的心理：你用了程序就要老老实实引文章！（这也是真的）

希望各位能在更多实践中使用它，折磨它，给它找茬，并且希望它对你预测反应活性位点的时候有用。
如果你发现问题，请务必把bug、改进意见什么的都砸给我，联系方式在manual里
论文刚被RSC的《React. Chem. & Eng.》接收（https://doi.org/10.1039/D5RE00506J），编辑部还没排版，各位凑合着看——其实我觉得也够清楚了——如果真的把排完版后的pdf传上来，怕是还要遇到版权问题呢。本帖的附件里就是论文本身（算法、推导过程、实例测试、软件源码、软件说明书）
其实，一个纯的算法论文，能发在if=3点几的杂志上，也行，不亏不赚吧
这个程序其实来源于我2019年的一个脑洞（http://bbs.keinsci.com/thread-11899-1-1.html）。后来断断续续想了一段时间，毕竟我的主业也不是搞算法。现在终于落地了。
总之，请大家试用，找茬，我改进，一条龙。

下面是一些示例（被一层球壳儿包着的，就是反应活性位点。球壳儿本身就是被进攻的“来源方向”）

图1：Ru被进攻的方向，残缺的壳儿表达了进攻方向以外侧为主，但也不排除可以从内侧的缝儿里挤进去

Ru的被进攻的来向，主要是顶头儿的外侧

图2：B被进攻的方向，上下两个壳儿，表达了进攻方向是在BR3平面的上下两侧。

B的被进攻的来向，是BR3的正上方和正下方

图3：360度环形进攻中间的C，但不进攻左边的C或者右边的O

中间的C，环绕进攻都行

图4：对O而言，这是非常清晰的外侧进攻方式

O的被进攻的来向，是外侧

图5：极其清晰的表明了卡宾C的进攻方向

卡宾的C的被进攻的来向

”

图6：竞争性反应的预测也不在话下


图7：真不怕双位点。其实这种靠得很近的O和N，如果来了个氢离子的话，氢离子本来就是在O和N之间来回跳转着玩儿，根本稳定不下来。



图8：这个其实很有意思。你把isovalue设为2亿的时候，你发现N似乎是被360度环形进攻的。但你把isovalue设为200亿的时候，你会发现优势进攻位点其实是分子平面的两侧。
这个反应还有个插曲。审稿人问，为啥N的孤对电子处却被不能被进攻？我当时心里吐槽：大哥你疯了吧？这是亲核进攻，不是亲电进攻，当然要躲着N的大大的孤对电子啦！


各位不用怕我“挑反应”。这个算法在验证阶段的时候，我和研究生是直接拿了一本纸质的J. J. Li写的《Name Reactions: A Collection of Detailed Mechanisms and Synthetic Applications》，我们随机抽签儿，抽到哪个算哪个，就是为了保证随机性。

不过呢，最大的问题是我们一共验证了二十几个反应而已。到你们手里，怕是能验证出一大堆问题来。啊不，我不怕你们验证出问题来，我期待你们验证出问题来。谁验证出问题，谁联系我，将来出version2的时候，我把你们列到acknowledge或者author list里。不过说好了啊，只提交bug，只能进acknowledge；一起合作改进，才能进author list。

在version2里，我要改的东西包括：
(1) 我自己发觉的：现在算出来的数值呢，是个正数，但是这是一个从极小到极大的正数，小可能是十的负十几次方，大可能是10的几百次方。所以我在version2里打算用Sigmoid类归一化：1/y=1+exp[c-kln(x)]。暂时来说，c=0，k=0.01吧。
(2) 也是我自己发觉的：现在的原子实，其实是用原子半径乘以缩放系数来模仿的，生硬了。我觉得下一版本，应该改成真正的，用波函数考量的真正原子实。
(3) 还是我自己发觉的：现在的version1其实是锁死了298K。version2里，我得考虑变温。不过现在也能用。
(4) 审稿人发掘的：软硬酸和软硬碱，在这里能考虑吗？怎么考虑？这是个大问题，但我没思路。但我承诺审稿人在version2里解决。（啊，疯！）
(5) 大家继续……

wal

程序直接附在论文里好啊 赞
见过发论文挂GitHub仓库假装开源的那种 发完论文仓库直接撤了 令人感到幽默

yjcmwgk

wal 发表于 2026-2-20 18:51
程序直接附在论文里好啊 赞
见过发论文挂GitHub仓库假装开源的那种 发完论文仓库直接撤了 令人感到幽默

哈哈哈哈哈哈哈哈 这是我的习惯
如果论文里我用的什么自编代码 我总是随着论文一起发表
明面上的理由：为了让大家复现计算结果更方便！
暗戳戳的心理：你用了我的程序就要老老实实引我文章！

牧生

由于我们没有买gaussian，所以我用了orca来计算了单点能，以期望重复楼主的操作。以我的理解，主要就是要得到某个分子0电荷单点能和静电势，+1和-1电荷激发态的单点能。

但有一步卡住了，请楼主帮忙解答一下。

为了在我的笔记本电脑上操作，使用的计算级别比较低。我用EthyleneKetene进行计算

1. 在b97-3c级别下，优化电中性的EthyleneKetene，并得到了当前级别下的单点能文件EthyleneKetene.gbw，直接用Multiwfn自带的ESPiso.bat （计算ESP时改成了高质量格点的选项）计算得到了ESP1.cub和density1.cub文件，分别改名为ESP0.cub和density0.cub

2. 用优化后的EthyleneKetene，再次生成单点能计算文件，电荷和自旋多重度该为1 2，同样方法得到density+1.cub

3. 改2中的电荷和自旋多重度该为-1 2，同样方法得到density-1.cub

4.把上述得到的4个cub文件和RAI.exe放在同一个文件夹里面，打开RAI.exe, 依次把density-1.cub，density0.cub，density+1.cub，ESP0.cub拖进去    （拖的顺序是带负电，0电荷，正电荷）

---- RAI(+) and RAI(-) Calculating Program (v 1.0) ----
Input the cube file name of negative charge file ( -1 e ) unit of length must be Bohr, the same below:
C:\Users\gxw85\Desktop\g6\1\density-1.cub
Input the cube file name of electron neutral file (  0 e ):
C:\Users\gxw85\Desktop\g6\1\density0.cub
Input the cube file name of positive charge file  ( +1 e ):
C:\Users\gxw85\Desktop\g6\1\density+1.cub
Input the cube file name of electrostatic file    (  esp ):
C:\Users\gxw85\Desktop\g6\1\ESP0.cub
ATTENTION!!! Cubes definitions are different, please check!
PAUSE
To resume execution, type go.  Other input will terminate the job.



在这里就不对了。请问可能错在哪一步呢？

yjcmwgk

牧生 发表于 2026-2-24 14:22
由于我们没有买gaussian，所以我用了orca来计算了单点能，以期望重复楼主的操作。以我的理解，主要就是要得 ...

就是这一句啊
ATTENTION!!! Cubes definitions are different, please check!

你去查cube文件的头信息，是不是不一样？

你自己去看看论文里的scheme 1
你卡在这儿，第一个判断框里
（我觉得自己好贴心啊，在论文里直接贴程序流程框图）

牧生

yjcmwgk 发表于 2026-2-24 15:38
就是这一句啊
ATTENTION!!! Cubes definitions are different, please check!

谢谢贴心指导，确实是因为这个原因。

现在已经完全走通了流程

yjcmwgk

牧生 发表于 2026-2-24 17:02
谢谢贴心指导，确实是因为这个原因。

现在已经完全走通了流程

啊 谢谢你的反馈 正好借楼重申一下 如果cube的头定义（格子定义）不一样 那两个cube互相运算将会毫无意义 因为每个空间采样点的坐标都不一样 完全没意义了 所以这个程序一开始就会检测cube的头定义