树枝状大分子的动力学模拟问题-结构优化 RESP sdf

wuyufei

自学了一段时间MD，目前合成了一系列含自由基（TEMPO）的树枝状大分子，想研究其EPR(Electron paramagnetic resonance)性质，想通过动力学模拟查看构型变化，水分子分布，帮助解释结构对epr曲线影响。典型的分子结构如下，末端含有TEMO自由基，中间链上为钠盐，水溶性。

原子数为    G0(314)      
              
原子数为 G3(2514)
以下是我的步骤 （314个原子G0）,使用GAFF力场，GROMACS做MD

• Chemdraw画出分子，适当结构透视，复制到Chem3d，MM2能量最小化。
• 使用XTB的GFN-FF力场优化。 xtb G0.xyz --gfnff --gbsa h2o --charge -6 --uhf 6 --opt
  

（问题2.1：尝试使用高斯进行结构优化，得到的结构看起来还行，但无奈很难能量收敛，对于再大的分子G3，SCF不收敛，不能进行opt。
# opt upm7

G0

-6 13
xyz
  问题2.2  使用MOPAC不会正确设置电荷和spin， 得到的结构很怪，比如苯环的氢不共面，使用PM6-D3H4得到结构错误更多。所以最终选择了XTB。
PM7 MOZYME CHARGE=-6 EPS=78.4
  问题2.3 关于XTB 6.5.1,对于G0（314个原子）可以正常opt，对于更多的原子如超过500，但更多时候会出现错误，不知道什么问题？
* calculating model hessian...
Using GFN-FF Lindh-Hessian
* fragmented diagonalization... 2 fragments
forrtl: severe (170): Program Exception - stack overflow
Image              PC                Routine            Line        Source
xtb.exe            00007FF742C2CF97  Unknown               Unknown  Unknown
xtb.exe            00007FF741076EA0  XTB_GFNFF_FRAGHES         328  frag_hess.f90
xtb.exe            00007FF740D5C6DD  XTB_RELAXATION_EN         645  relaxation_engine.f90
xtb.exe            00007FF740AD0F1C  XTB_GEOOPT_mp_GEO         102  geoopt_driver.f90
xtb.exe            00007FF74098892F  XTB_PROG_MAIN_mp_         643  main.F90
xtb.exe            00007FF7409971C7  MAIN__                     55  primary.f90
xtb.exe            00007FF742AB051E  Unknown               Unknown  Unknown
xtb.exe            00007FF742C2D1A0  Unknown               Unknown  Unknown
KERNEL32.DLL       00007FFB6F6F7034  Unknown               Unknown  Unknown
ntdll.dll          00007FFB70FE26A1  Unknown               Unknown  Unknown


3.使用g16计算RESP。
# ub3lyp/6-31g(d) scrf=(smd,solvent=water) pop=mk empiricaldispersion=gd3 iop(6/33=2,6/42=6,6/50=1)
  问题3.1 RESP应该对结构影响很大吧，直接计算从XTB得到的结构会不会有问题？上述关键词合理吗？
  问题3.2 无法使用常见的在线工具获得拓扑文件，怀疑应该是自由基的问题。

4. 使用AmberTools加载电荷、力场、电荷和水模型。
antechamber -i G0.qfi -fi gout -o G0.mol2 -fo mol2 -pf y -c resp
parmchk2 -i G0.mol2 -f mol2 -o G0.frcmod
tleap
source oldff/leaprc.ff14SB
source leaprc.gaff
loadamberparams G0.frcmod
G0=loadmol2 G0.mol2
source leaprc.water.tip3p
solvatebox G0 TIP3PBOX 10
addions G0 Na+ 0
check G0
quit
问题4 这里使用GAFF是否合理？
问题5 如果计算在二氯甲烷中的MD，加载DCM的拓扑文件，solvatebox G0 DCM 10。试了一下，dcm排列很规则，进行npt会出错，需要先Berendsen，再Parrinello-Rahman。
最后也能顺利MD。这样得到的结果有没有意义？或者直接加载优化好的溶剂盒子，但还没有学会怎么弄拓扑文件。

5. acpype将拓扑文件转换为gromacs格式。 acpype -p G0.prmtop -x G0.inpcrd -d
6. 使用GROMACS 2019.2 GPU进行MD。能量最小化，NVT 1ns, NPT 1ns,  MD 100ns。温度密度压力等均稳定。
7. 数据分析，
居中，消除平动转动，gmx spatial计算了10ns MD 中心分子周围的水分子分布.

vmd查看，isovalue0.02，中间水分子较多
  问题7 不同分子的sdf的是否可以比较，比如相同的isovalue，A分子的水看起来比B的多，那么能说A的周围出现的水更多吗？

   这样一套下来，流程是否合理，得到的结果可信吗（能发文章吗）？对于此任务老师们有没有别的看法
谢谢各位老师，以上均为查论坛资料加自己领悟的，肯定有很多不足，想听听大家关于上述问题的看法，谢谢啦~

wuyufei