帮助设定Gaussian输入文件中优化flag和MM电荷的小工具

sobereva

帮助设定Gaussian输入文件中优化flag和MM电荷的小工具

A small tool to facilitate setting up optimization flags and MM charges in Gaussian input file

文/Sobereva @北京科音   2016-Jun-6

近来由于计算需要写了两个小程序setopt和setcharge，一个是用来帮助设定Gaussian输入文件中的优化flag，另一个是帮助设定MM电荷。虽然程序很简单，但很有用，这里分享下。Windows版及源代码下载地址： optflag 1.0.rar (336.18 KB, 下载次数 Times of downloads: 78)

2016-6-6 12:19 上传 Uploaded

点击下载
Click to download

setcharge 1.0.rar (415.66 KB, 下载次数 Times of downloads: 83)

2016-6-6 12:19 上传 Uploaded

点击下载
Click to download

1 optflag程序

Gaussian里optimization flag就是原子后面写0、-1来分别代表优化时允许移动和被冻结。对于ONIOM，设成同样的负值的原子还可以被当成刚性片段优化。optflag小程序可以将指定的原子的优化flag设为指定值。来看一个例子:



这是个靛蓝(indigo)的分子团簇，从晶体结构中截取出来的。pdb文件在此： cluster.pdb (353.98 KB, 下载次数 Times of downloads: 47)

2016-6-6 12:18 上传 Uploaded

点击下载
Click to download

我们想把中间红色的分子的周围的一圈分子，即黄色区域的靛蓝分子的优化flag设为0，允许它们在优化时调整结构，而将最外围那些以及中间的红色靛蓝分子设为冻结。

要做到这个目的，我们先用VMD生成原子序号列表。选择范围对应于：
中心分子：residue 43
中心分子+相邻部分：same fragment as within 6 of residue 43

比如要生成中心分子的序号列表，就在VMD控制台里输入atomselect top "residue 43"，比如提示atomselect0，就再输入atomselect0 list得到序号列表和atomselect0 num得到原子数。然后写成索引文件center.txt，将用于作为optflag的输入，内容如下。
-30
1290 1291 1292 1293 1294 1295 1296 1297 1298 1299 1300 1301 1302 1303 1304 1305 1306 1307 1308 1309 1310 1311 1312 1313 1314 1315 1316 1317 1318 1319

第一行是原子数，之后是原子序号。由于VMD输出的序号是从0开始的，原子数要写成负值来指明这一点。以同样步骤，生成中心+相邻部分分子的索引文件small.txt。

然后用gview打开.pdb文件，保存成.gjf文件，把原子坐标那部分截出来，并且去掉(pdbname...)那一堆，保存到all.txt里面。

双击启动optflag，输入
all.txt
4530   //总原子数
0   //当前输入文件里还没有优化flag
all   //选全部原子
-1   //要设成的优化flag
此时提示在当前目录产生了new.txt，打开一看，已经相对于all.txt把-1的优化flag加上了。

然后把中心+相邻部分分子优化flag设为0，启动optflag输入
new.txt
4530
1   //输入文件里已经有优化flag了
small.txt   //用来选定中心分子+相邻部分的索引文件
0   //要设的优化flag

最后一次，将中心分子优化flag设为-1，启动optflag输入
new.txt
4530
1
center.txt
-1

好啦，此时得到的new.txt的内容就已经可以拷到Gaussian输入文件的坐标部分了，优化时就只有中心分子临近的一圈分子会被优化了。感兴趣的话可以把gjf文件拉进gview，进edit-atom groups，选freeze分类，然后高亮或隐藏不同的部分检查是否设定合理。


2 setcharge程序

Gaussian中用分子力学计算的时候涉及到给原子设定原子电荷，对于体系中包含许多同种分子，比如包含一大堆水分子，或者上面的例子包含一堆indigo，要给自行挨个给同种分子填上原子电荷信息很麻烦。setcharge可以帮助做这个事情。

对于上面的例子，我们要给团簇所有indigo分子设定在indigo孤立状态下算出的CHELPG电荷，并且假定中间的indigo分子带+1电荷，因此要给它赋上indigo阳离子状态下计算的CHELPG电荷。做法如下：

先对indigo在中性和+1状态下用Multiwfn计算CHELPG原子电荷（参考http://sobereva.com/441），然后分别写到indigo.txt和indigo+1.txt，第一行是原子数，之后是每个原子的电荷，例如indigo.txt
30
0.099919
0.519588
-0.183248
-0.016362
-0.189175
0.020601
-0.308034
0.418772
-0.659064
-0.562279
0.099919
0.519588
-0.659064
-0.183248
-0.562279
0.418772
-0.016362
-0.308034
-0.189175
0.020601
0.102237
0.100428
0.081598
0.143281
0.431738
0.102237
0.100428
0.081598
0.143281
0.431738

先给所有indigo都赋上中性下indigo的CHELPG电荷，启动setcharge，输入
new.txt  //这个是上个例子最后得到的new.txt
4530
1
all   //给全部原子赋上电荷
indigo.txt  //含有原子电荷的文件

原先new.txt里的优化flag还都保留着。当前体系一共4530个原子，indigo.txt里有30个原子的电荷，所以会给1~30、31~60、61~90...4501~4530都依次赋上indigo.txt里的电荷。此时可以看到new.txt里的原子变成了诸如这样
C---0.308034
即曰C的原子电荷为-0.308034。前两个横杠之间是原子类型，这里没设所以留空。

然后给中心原子赋上+1态的电荷。启动setcharge，输入
new.txt  //这个是上个例子最后得到的new.txt
4530
1
center.txt  //中心分子的索引文件，见上例
indigo+1.txt  //含有+1态CHELPG电荷的文件

然后new.txt里的信息就可以拷到gjf文件里用了。

小范范1989

sob老师好，通过老师这个帖子学到了很多，对我之前的一个想法帮助非常大，（想法之前和您讨论过：ONIOM优化，但是在算频率的时候，把MM区域当做背景电荷处理）

所以，针对您的第二个方案setcharge，我这里有几个问题想咨询一下老师：
1：我看老师说这个CHELPG电荷，会给后面晶体中每个分子这个附上电荷。我想咨询的是高斯计算的CHELPG电荷的顺序，和咱们要对晶体中赋值的顺序是一样的吗？如何才能保证一样？
其实我的意思是，会不会出现高斯计算的H的电荷被赋值到了晶体中C原子上面。这个对应关系如何保证准确？
2：希望sob老师指点一下我计算频率，把MM当做背景电荷的处理方式：
  2.1：通过老师的optflag程序，我对晶体结构进行了固定等等，然后用ONIOM进行了优化，得到了优化的坐标文件。
  2.2：现在，我想通过老师说的setcharge方式来把我周围的分子做背景电荷，
         处理方式为：先通过老师的这个setcharge程序对所有的分子进行电荷赋值，
         然后：把中心分子的坐标单独取出来，保存做高斯输入文件（去除刚刚计算的电荷）。同时把周围所有分子的坐标和刚刚得到的电荷，做背景电荷处理，也就是放到中心分子坐标的后面。


3：不知道我2.2这种方式对不对，也就是中间分子做中性分子，不考虑自身的电荷，而把周围的做背景电荷。但是我看sob老师是把中间的做成了带+1的文件，这里老师只是为了展示计算还是说ONIOM处理的时候的一种方式？谢谢老师。

sobereva

小范范1989 发表于 2016-6-7 08:16
sob老师好，通过老师这个帖子学到了很多，对我之前的一个想法帮助非常大，（想法之前和您讨论过：ONIOM优化 ...

1 团簇里相同种类分子原子顺序通常都是一样的，你取其中一个，优化，算电荷，显然给出来的顺序还是和团簇里一致

2.2 你没必要去除中心分子的电荷，你不用分子力场来做的话原子电荷设定就不会被利用而已，不碍事。

3 对。+1只是为了演示，另外涉及到算外重组能的时候也需要这么干。

小范范1989

sobereva 发表于 2016-6-7 09:19
1 团簇里相同种类分子原子顺序通常都是一样的，你取其中一个，优化，算电荷，显然给出来的顺序还是和团簇 ...

非常感谢sob老师，学会了好多。
PS：顺带一下，这个gromacs培训班是暑假开学之后（9月后）？谢谢老师

sobereva

小范范1989 发表于 2016-6-7 09:35
非常感谢sob老师，学会了好多。
PS：顺带一下，这个gromacs培训班是暑假开学之后（9月后）？谢谢老师

是

小范范1989

sobereva 发表于 2016-6-7 09:19
1 团簇里相同种类分子原子顺序通常都是一样的，你取其中一个，优化，算电荷，显然给出来的顺序还是和团簇 ...

再咨询sob老师一下，这里为什么使用CHELPG电荷呢？原因是什么，或者参考文献，因为我怕后面写的时候审稿人问，就想先说一下为什么用CHELPG电荷。谢谢老师第二个就是：计算的高斯输出文件，读取电荷是哪一部分呢？我看有         
Condensed to atoms (all electrons):
Mulliken charges:
            .......
Mulliken charges with hydrogens summed into heavy atoms:


Fitting point charges to electrostatic potential
Charges from ESP fit, RMS=   0.00431 RRMS=   0.38174:
ESP charges:


ESP charges with hydrogens summed into heavy atoms:

以及
-----------------------------------------------------------------

              Electrostatic Properties (Atomic Units)

-----------------------------------------------------------------
    Center     Electric         -------- Electric Field --------
               Potential          X             Y             Z
-----------------------------------------------------------------
  这么多选项呢，

PS：刚刚qq群里，sob老师已经指点，是读取Charges from ESP fit, RMS=   0.00431 RRMS=   0.38174:
ESP charges:
谢谢sob老师

小范范1989

sobereva 发表于 2016-6-7 09:19
1 团簇里相同种类分子原子顺序通常都是一样的，你取其中一个，优化，算电荷，显然给出来的顺序还是和团簇 ...

sob老师好，现在遇到这个背景电荷赋值的时候和初始分子坐标不对应的问题。
我的分子如附件。
主要是我分子的编号问题，同一个分子内，虽然编号连续了，但是相同分子之间，编号并不是相同的，比如说第一个是S-N-C-C-H，而第二个分子却是S-C-N-H-C这个样子，不知道如何把这些分子的顺序一致的排列起来。
希望sob老师能帮助调节一下这两个分子，能够指点一下如何实现，谢谢sob老师。
PS：之前用sob老师的这个方式完全没有问题。我感觉这次是因为实验分子给的这个cif的原因、cif少了H，自己扩充了层数，然后通过又保存坐标，同时添加了H。所以出现了顺序不对应的问题。
谢谢sob老师

sobereva

小范范1989 发表于 2016-6-30 17:00
sob老师好，现在遇到这个背景电荷赋值的时候和初始分子坐标不对应的问题。
我的分子如附件。
主要是我 ...

这个得自己写个程序，否则没什么办法

小范范1989

sobereva 发表于 2016-7-1 01:31
这个得自己写个程序，否则没什么办法

偶，这样子，好的， 谢谢sob老师。

rendong

老师您好，我想请教一下，在Gaussian中，使用ONIOM对有机分子晶体进行激发态结构优化，低层部分使用PM7，是否需要给低层部分进行电荷赋值

sobereva

rendong 发表于 2021-7-7 15:40
老师您好，我想请教一下，在Gaussian中，使用ONIOM对有机分子晶体进行激发态结构优化，低层部分使用PM7，是 ...

不需要
只有力场才需要原子电荷

rendong

sobereva 发表于 2021-7-8 12:27
不需要
只有力场才需要原子电荷

谢谢老师，这下我就不用重新再跑一次了