如何使用Sobtop创建铁硫蛋白中铁硫簇[Fe2S2]的GROMACS的拓扑文件

MaYong

各位老师好，

我最近在进行含一个 [Fe₂S₂] 铁硫簇蛋白的分子动力学（GROMACS）模拟，在铁硫簇参数化过程中遇到比较棘手的问题，特来请教。

主要参考了以下两种方法：（1） http://sobereva.com/soft/Sobtop/#FAQ 中FAQ7中的操作步骤：“比如蛋白质中有个铁硫簇和几个残基结合，那就把铁硫簇以及与之配位的残基抠出来并饱和边缘作为一个簇模型，冻结边缘原子，用量子化学程序做优化和振动分析”。（2）http://bbs.keinsci.com/thread-27803-1-1.html中的使用Sobtop超级方便地创建二茂铁的GROMACS的拓扑文件  这两种方法进行操作。
碰到一些问题：
我首先提取了FE2S2以及与两个FE配位的4个CYS作为一个簇,然后对CYS上的除了S原子以外的所有原子进行了加氢处理以使其饱和。然后在跑gaussian的时候用了教程中的计算级别TPSSh/def2-TZVP。但是过程中SCF始终不收敛，无法达到收敛标准，多次尝试解决后问题仍然存在。

sobereva

上传Gaussian输入输出文件

MaYong

sobereva 发表于 2025-12-18 05:33
上传Gaussian输入输出文件

gaussian 输入文件：

%mem=64GB
%chk=F2R_4CYS.chk
%nprocshared=72
#p opt freq def2tzvp tpssh

F2R_4CYS

-3 2
N                  2.67658132   -3.27384048   -1.37845297
H                  2.44687145   -3.73010534   -0.51876945
C                  1.83036363   -3.80098435   -2.45257212
H                  2.33973652   -4.66703175   -2.87783106
C                  0.50616418   -4.29811432   -1.86434036
H                 -0.04658333   -4.84481481   -2.63005751
H                  0.71119442   -4.97175448   -1.03118879
S                 -0.51195912   -2.90326066   -1.27976426
C                  1.57131392   -2.81471454   -3.60109427
O                  1.31780483   -3.26956166   -4.71098626
N                 -2.72747278   -3.43702568    1.50152288
H                 -2.28194725   -2.77318154    0.90084164
C                 -1.95104924   -3.57759045    2.73426090
H                 -2.50018177   -4.22578978    3.41734150
C                 -0.59702831   -4.25380182    2.42916832
H                 -0.20295680   -4.68992219    3.34809142
H                 -0.74754847   -5.05976683    1.70939747
S                  0.63547336   -3.08173622    1.77956544
C                 -1.77534515   -2.22427318    3.43042499
O                 -2.05483054   -1.18505349    2.84799444
N                  1.94738015    1.67974445    3.41944115
H                  1.60477196    1.14801475    2.64491991
C                  1.94541260    3.10167378    3.07389070
H                  2.92935946    3.51260937    3.27879682
C                  1.68360169    3.23482811    1.56355598
H                  1.85880573    4.25617166    1.24297074
H                  2.37613884    2.58952862    1.02221426
S                 -0.01325240    2.74342333    1.17748146
C                  0.96823012    3.94475389    3.91409505
O                  0.57889325    5.03040493    3.48661218
Fe                 0.07380736    0.07173835    1.48869353
Fe                -0.25426158   -0.31745101   -1.22653453
S                  1.66619525   -0.10328617   -0.13797952
S                 -1.89508278   -0.03848664    0.36764666
N                 -0.56967783    5.05256612   -1.50499041
H                 -0.61385223    4.51198319   -0.66486023
C                 -0.41265034    4.38409335   -2.79837029
H                  0.27525617    4.96809099   -3.40270495
C                  0.22899914    3.00030386   -2.58813570
H                  0.59105195    2.63838005   -3.55107595
H                  1.07139341    3.08868957   -1.90587705
S                 -0.92696730    1.77756266   -1.90455719
C                 -1.71371413    4.25285397   -3.62212041
O                 -1.65276258    3.77531873   -4.75188100
H                 -0.63273641    6.04964442   -1.46188050
H                  2.87988247    1.39058727    3.63583788
H                  3.63755961   -3.44156082   -1.59843163
H                 -2.78122236   -4.32185267    1.03871387


gaussian 输出文件：

Cycle 128  Pass 1  IDiag  1:
RMSU=  2.22D-07    CP:  1.00D+00  9.77D-01 -1.04D-02  1.70D+00  2.39D+00
                    CP: -8.79D-02  1.01D+00
E= -5906.15722508253     Delta-E=        0.000326100233 Rises=F Damp=F
DIIS: error= 1.98D-03 at cycle 128 NSaved=  20.
NSaved=20 IEnMin= 3 EnMin= -5906.16984338742     IErMin= 1 ErrMin= 5.08D-04
ErrMax= 1.98D-03  0.00D+00 EMaxC= 1.00D-01 BMatC= 2.04D-02 BMatP= 3.41D-03
IDIUse=3 WtCom= 9.95D-01 WtEn= 5.08D-03
Rare condition: small coef for last iteration:  0.000D+00
Coeff-Com: -0.261D+00 0.110D+00-0.606D+00 0.923D-01 0.731D+00 0.627D-01
Coeff-Com: -0.503D+00-0.975D+00 0.136D+01 0.218D+01-0.185D+01 0.539D+00
Coeff-Com:  0.537D+00-0.400D-01 0.260D+00-0.444D+00 0.343D+00 0.407D+00
Coeff-Com:  0.556D+00-0.150D+01
Coeff-En:   0.000D+00 0.000D+00 0.000D+00 0.000D+00 0.000D+00 0.000D+00
Coeff-En:   0.000D+00 0.000D+00 0.000D+00 0.000D+00 0.000D+00 0.000D+00
Coeff-En:   0.000D+00 0.000D+00 0.000D+00 0.000D+00 0.000D+00 0.000D+00
Coeff-En:   0.000D+00 0.100D+01
Coeff:     -0.261D+00 0.110D+00-0.606D+00 0.923D-01 0.731D+00 0.627D-01
Coeff:     -0.503D+00-0.975D+00 0.136D+01 0.218D+01-0.185D+01 0.539D+00
Coeff:      0.537D+00-0.400D-01 0.260D+00-0.444D+00 0.343D+00 0.407D+00
Coeff:      0.556D+00-0.150D+01
Gap=     0.015 Goal=   None    Shift=    0.000
Gap=     0.014 Goal=   None    Shift=    0.000
RMSDP=6.52D-05 MaxDP=1.01D-02 DE= 3.26D-04 OVMax= 4.17D-02

>>>>>>>>>> Convergence criterion not met.
SCF Done:  E(UTPSSh) =  -5906.15722508     A.U. after  129 cycles
            NFock=128  Conv=0.65D-04     -V/T= 2.0022
<Sx>= 0.0000 <Sy>= 0.0000 <Sz>= 0.5000 <S**2>= 5.6670 S= 1.9325
<L.S>=  0.00000000000   
KE= 5.893267313231D+03 PE=-2.470285755063D+04 EE= 7.607241670462D+03
Annihilation of the first spin contaminant:
S**2 before annihilation     5.6670,   after    17.8652
Convergence failure -- run terminated.
Error termination via Lnk1e in /public1/home/sch3181/software/g16/l502.exe at Thu Dec 18 15:36:21 2025.
Job cpu time:       3 days 14 hours  7 minutes 42.3 seconds.
Elapsed time:       0 days  5 hours 26 minutes 36.2 seconds.
File lengths (MBytes):  RWF=   1023 Int=      0 D2E=      0 Chk=     35 Scr=      1

sobereva

MaYong 发表于 2025-12-18 15:49
gaussian 输入文件：

%mem=64GB

当前结构完全离谱，截断的地方都没加氢饱和，SCF容易收敛就怪了
优化目的Fe用SDD、配体6-31G*足矣。def2-TZVP单点计算才用得着
当前体系要当反铁磁性耦合体系算，看
谈谈片段组合波函数与自旋极化单重态
http://sobereva.com/82