Study note：构象搜索

在实验中观察到了奇葩现象，需要通过计算验证形成exciplex的可能性。由于此前半年的学习中一直专注于单分子的计算，在上手exciplex时建模遇到了一些麻烦，作本文聊以记录备忘。

构象搜索

对于团簇或柔性体系，建模是个大麻烦。手动摆的位置很可能优化到局域极小点而非全局极小点，此时就需要molclus结合这些工具来进行构象搜索：

• 团簇构象搜索：genmer。
• 分子构象搜索：gentor不支持环状区域构象搜索，Confab或Frog2更简单但不支持环状区域构象搜索，不支持有机分子以外的情况。
• md采样：既可用于团簇也可用于分子。gromacs等使用繁琐，需要有合适的力场；xtb做md效率高，普适。共同缺点可能有采样不足的问题。

genmer

genmer是molclus工具包中提供的生成团簇初始构象的工具。以二聚体初始构象为例，settings.ini设置如下：

  ngen= 50   // The number of configurations expected to be generated. Default: 1
  ngenmust= 50  // The number of configurations must to be generated. Default: 0 (No requirement)
  step= 0.1  // Step size to grow monomer from cluster center, larger value leads to faster   calculation but lower accuracy. Default: 0.1
  ishuffle= 0    // 0: Don't shuffle the sequence of monomer types  1: Shuffle the sequence. Default: 0
  imoveclust= 1    // 1: Move geometry center of the constructed clusters to (0,0,0)  0: Keep initial position. Default: 1
  irandom= 1    // 0: Result of each time of execution will be identical  1: Will be different to due use of different random seed in each run. Default: 1
  iradtype= 1    // 0: Use CSD covalent radii to detect contacts  1: Use Bondi vdW radii. Default: 1
  radscale= 1.0  // Scale factor of atomic covalent and vdW radii. Default: 1.0
  maxemit= 300  // The maximum number of attempts of emitting monomer. Default: 300
---- Below are number (may with options) and .xyz file of each type of monomer ----
2 noxrot noyrot rmax=5
bodipy-dimer/IPR.xyz

由于绝大部分地方采用了默认设置，可以简化为：

ngenmust= 50
----
2 noxrot noyrot rmax=5
dimer/bodipy.xyz

由于genmer运行极快，在没有私人计算节点的情况下，完全可以用个人电脑运行genmer来生成初始构象。将bodipy.xyz放在genmer.exe目录下，运行完成后生成了traj.xyz，可以在VMD中观看生成的结构。

运行后，使用isostat.exe进行聚类，最终得到cluster.xyz文件。

Ref：

genmer：生成团簇初始构型结合molclus做团簇结构搜索的超便捷工具。

molclus

genmer生成的初始构象可以用molclus来处理。molclus的配置文件为：

  iprog= 4  // The computational code to invoke. 1: Gaussian, 2: MOPAC, 3: ORCA, 4: xtb, 5: Open Babel
  ngeom= 0  // Geometries in traj.xyz to be considered. 0: All, n: First n. You can also use e.g. 2,5-8,12-14 to customize range
  itask= 0  // Type of task. 0: Optimization, 1: Single point energy, 2: Frequency, 3: Optimization with frequency, -1: Composite method (only for Gaussian, e.g. CBS-QB3)
  ibkout= 0  // When backup output file. 0: Never, 1: All, 2: Successful tasks, 3: Failed tasks, -1: If backup file exists, directly loading  rather than calculating (only available for Gaussian and ORCA)
  distmax= 999  // If distance between any two atoms is larger than this (Angstrom), the geometry will be skipped. If =0, load maximum distance between atomic pairs from distmax.txt in current folder
  ipause= 0  // When pause molclus. 1: Optimization didn't converge (For Gaussian and ORCA), 2: After each cycle, 0: Never pause
  iappend= 0  // 0: New systems will replace the old isomers.xyz, 1: New systems will be appended to the old isomers.xyz
  freeze= 0  // Index of atoms to be freezed during optimization, e.g. 2,4-8,13-14. If no atom is to be freezed, this parameter should be set to 0
--- Below for Gaussian ---
  gaussian_path= "g16"  // Command for invoking Gaussian
  igaucontinue= 1  // 1: If optimization exceeded the upper number of steps, continue to optimize the last geometry using template2.gjf and then template3.gjf 0: Don't continue
  energyterm= "HF=" // For itask= 0 or 1, set label for extracting electron energy from archive part of output file, e.g. HF=, MP2=, CCSD(T)=; if set to "TD=", will load TDDFT excited state energy of interested state. For itask= -1, this sets label for extracting thermodynamic quantity, e.g. "G4 Free Energy=" means G4 free energy will be extracted (need to write double quotes!). For itask= 2 or 3, this parameter does not need to be set because free energy is always extracted (however, if template_SP.gjf is presented, this parameter also determines loading which energy from high-level single point task). 
  ibkchk= 0  // The same as ibkout, but for .chk file
--- Below for ORCA ---
  orca_path= "/sob/orca/orca"  // Command for invoking ORCA
  ibkgbw= 0  // The same as ibkout, but for .gbw file
  ibktrj= 0  // 1: Backup optimization trajectory (*_trj.xyz) with step number in the file name, 0: Don't backup
  ibkhess= 0  // 1: Backup .hess file with step number in the file name, 0: Don't backup
  mpioption= none  // Special MPI option for running ORCA. e.g. --allow-run-as-root. If no option is needed, this option should be set to "none"
--- Below for MOPAC ---
  mopac_path= "MOPAC2016.exe"  // Command for invoking MOPAC
--- Below for xtb ---
  xtb_arg= "--gfn 2 --chrg 0 --uhf 0"  // Additional arguments for xtb, e.g. "vtight --gfn 1 --gbsa h2o --chrg 0" (don't write task keywords such as "--opt" here)
--- Below for Open Babel ---
  obabel_FF= MMFF94  // Force field employed by Open Babel, available: MMFF94, Ghemical, UFF, GAFF
  obabel_param= "--steps 2500" // Additional arguments for Open Babel. Hint: For large system, "--log" could be added to monitor minimization process
--- Below for Shermo ---
  Shermo_path= "D:\CM\my_program\Shermo\Shermo.exe" // Command of invoking Shermo for calculating thermodynamic data for Gaussian and ORCA
  T= 298.15  // Temperature (K)
  P= 1.0  // Pressure (atm)
  sclZPE= 1.0  // Frequency scale factor for ZPE
  ilowfreq= 2  //Treatment of low frequencies. 0: Harmonic approximation. 1: Raising frequencies lower than 100 cm^-1. 2: Grimme's entropy interpolation

如果有个人的计算节点，将traj.xyz复制到molclus目录下，配置好后直接运行

nohup ./molclus > out.txt &

即可开始计算。而如果你像笔者一样需要提交到集群，就要写提交脚本了。

xtb

xtb运行也奇快无比，甚至跑100个构象这种任务用个人电脑也不是不行，只能说Grimme真是奇才。

# xtb的配置文件
source $HOME/apprepo/xtb/v6.7.0-intelmpi2021/scripts/env.sh
# 将 xtb 所在目录加入 PATH
export PATH=$PATH:~/apprepo/xtb/v6.7.0-intelmpi2021/app/bin
molclus的路径
export PATH=$PATH:~/apprepo/molclus_1.12_Linux

gaussian

待补充

xyz2QC

如果你是用的是集群，molclus的串行式运行可能不能把资源充分利用起来。可以使用xyz2QC.exe将traj.xyz拆分成Gaussian输入文件，然后提交到集群。坏处是后续能量提取等需要自己手动进行。