摘要：Flare V6的一个新特性是引入了配体QM计算。本文以(S)-Dexmedetomidine为例介绍了如何通过Flare调用PSI4进行ECD光谱计算，具体步骤包括：XED力场构象优化、XED力场构象搜索、QM构象优化、构象去重、激发态计算以及ECD图谱绘制。

肖高铿/2022-09-22

简介

Flare V6是一款集成了基于结构、基于配体与QM计算的综合性软件包，并引入多个新的特性[1,2,3]，其中一个新特性是配体分子的QM计算。之前已经介绍过，如何在Flare里用QM研究配体的静电特性[1]与构象稳定性[4]。本文以(S)-Dexmedetomidine(图1化合物)为例介绍如何通过Flare的Jupyter界面里调用PSI4[5]进行ECD光谱计算，详细的PSI4计算ECD图谱还可参见PSICON 2020的TDDFT教程[6]。本文的目的是演示，而非科研，因此更详细的ECD光谱计算规范请参考Pescitelli等人的文献[7]。

用Flare QM计算化合物的ECD图谱-墨灵格的博客

图1. (S)-Dexmedetomidine

Dexmedetomidine的SMILES代码如下：

1

C[C@@H](c1cccc(C)c1C)C2=CNC=N2

XED力场优化与构象搜索

首先，我们用Flare | Ligand | New创建(S)-Dexmedetomidine的3D结构（图2 步骤1），并视觉检查确保立体化学信息是正确的。然后用Flare | 3D Pose | Minimize使用XED对创建的结构进行精确优化（图2 步骤2与3）。

用Flare QM计算化合物的ECD图谱-墨灵格的博客

图2. 创建3D结构与XED力场优化

有了XED力场优化过后的3D构象，我们开始用XED力场进行构象搜索。点击 Flare | 3D Pose | Conf Hunt & Align （图2 步骤4），弹出图3所示的构象搜索与叠合对话框，使用如下参数：

  • Conformation Hunt： Very accurate and slow
  • Alignment：No Calculation
  • Ligands to Conf hunt & align：Selected ligand

然后点击start开始构象搜索，见图3步骤5。

用Flare QM计算化合物的ECD图谱-墨灵格的博客

图3. Conf Hunt & Align对话框

构象搜索结束，用Flare | 3D Pose | Conf Explorer可以发现总共找到了12个低能构象（能量窗=3kcal/mol，见图4左），这12个构象的侧视图与顶视图如图4所示。

用Flare QM计算化合物的ECD图谱-墨灵格的博客

图4. 12个低能构象

用QM优化构象

接下来用QM进一步对每个构象进行结构优化。点击Flare | 3D Pose| QM | Minimize Conformations（图5步骤6），弹出QM Minimize Conformers对话框，设置参数如下（图5步骤7）。优化分两个步骤进行：优化与单点能计算。

用Flare QM计算化合物的ECD图谱-墨灵格的博客

图5. QM优化构象

其中优化步骤参数如下：

  • Method： DFT
  • DFT functional：PBE
  • Use dispersion correction：Yes
  • Basis Set：6-31G(d)
  • Tolerance：normal
  • Max minimization steps：500

其中单点能计算步骤参数如下：

  • Method： DFT
  • DFT functional：B3LYP
  • Use dispersion correction：Yes
  • Basis Set：6-311+G(d,p)
  • Solvent：Methanol
  • PCM Model：IEFPCM

其中单点能计算的能量将用于构象系综分析，计算玻尔兹曼权重。溶剂的选择与实验数据采用的溶剂一致，也与下一步激发计算ECD图谱的溶剂一致。

计算完毕，点击Flare | 3D Pose | Conf explorer（图6步骤8）可以浏览每个优化的构象以及对应的XED与QM构象能。

用Flare QM计算化合物的ECD图谱-墨灵格的博客

图6. 构象分析Conf explorer对话框

在用QM优化之后，如图7所示，我们发现QM能量排序与XED力场能量排序不同，而且还有几个化合物的QM能量都为0。这提示，可能有部分构象优化之后得到相同的构象。为了减少后续TDDFT激发态计算的计算量，有必要将相同的构象合并（去重复）。

用Flare QM计算化合物的ECD图谱-墨灵格的博客

图7. QM构象能

在这个例子中，构象5与3重复，构象12与11重复，因此选择5与12（图8步骤8），先点击Del（图8步骤9），再点击accept则删除构象5与12（图8步骤10）。

用Flare QM计算化合物的ECD图谱-墨灵格的博客

图8. 构象去重复

最后，我们得到7个不重复的构象，将这些优化好的构象导出，保存为sdf文件备用。这个算例里，我将7个独立不重复的构象分别保存单个的SDF文件，共7个文件，分别为：CONF_1.sdf, CONF_2.sdf...,CONF_7.sdf。

激发态计算：ECD图谱生成

新建一个Flare project，点击Python | Jupyter打开Jupyter notebook。现在以最低能的构象CONF_1.sdf为例说明ECD光谱的计算，其它构象的ECD光谱计算采用完全一样的方法。

首先切换工作目录到保存构象的目录：

1
2
3

from cresset import flare
import os
os.chdir('/public/gkxiao/work/flare_ecd')

从SDF文件里提取分子的坐标，以便构造ECD计算的分子描述：

1
2
3
4
5
6
7
8
9
10
11

from rdkit import Chem
mol = Chem.MolFromMolFile('CONF_01.sdf',removeHs=False)
n = mol.GetNumAtoms()
formalcharge = str(Chem.rdmolops.GetFormalCharge(mol))
for i in range(n):
    pos = mol.GetConformer().GetAtomPosition(i)
    element = mol.GetAtoms()[i].GetSymbol()
    x = str(round(pos.x,4))
    y = str(round(pos.y,4))
    z = str(round(pos.z,4))
    print(element+" "+x+" "+y+" "+z)

输出如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31

C -0.9189 3.0772 0.5403
C -1.7467 4.0246 1.1222
N -1.5394 5.1889 0.4001
C -0.6108 4.9054 -0.5698
N -0.2164 3.6429 -0.5151
C -0.7283 1.6284 0.9049
C -1.2646 1.3199 2.317
C -1.3523 0.665 -0.1126
C -2.4974 1.0673 -0.8235
C -3.1342 0.2034 -1.7173
C -2.6248 -1.0845 -1.9131
C -1.4857 -1.5177 -1.2171
C -0.8422 -0.6433 -0.3022
C 0.3692 -1.1625 0.441
C -0.9554 -2.9135 -1.445
H 0.3649 1.4676 0.9034
H -2.4391 3.9861 1.9596
H -1.9863 6.0863 0.5731
H -0.2616 5.6557 -1.2793
H -1.0578 0.273 2.5925
H -2.3595 1.4573 2.352
H -0.8032 1.9801 3.0717
H -2.8794 2.0821 -0.6748
H -4.0222 0.5351 -2.2656
H -3.1126 -1.7694 -2.6162
H 1.1553 -1.498 -0.2602
H 0.8218 -0.4154 1.1075
H 0.1109 -2.0414 1.0609
H 0.0834 -2.903 -1.8244
H -0.9401 -3.5072 -0.512
H -1.5738 -3.4564 -2.1779

接着在B3LYP/6-311+G(d,p)理论水平下进行激发态计算，并以IEFPCM为溶剂模型考虑甲醇的溶剂效应：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73

import psi4
import numpy as np
from psi4.driver.procrouting.response.scf_response import tdscf_excitations
 
# Output file path
# arg 0: Path of output file
# arg 1: Whether overwrite or not
psi4.core.set_output_file("CONF_1_tddft.out",True)
 
# set threads and memory
psi4.core.set_num_threads(24)
psi4.set_memory(int(24e9))
 
# Scratch directory path
psi4.core.IOManager.shared_object().set_default_path("./tmp")
 
 
pcm_string = """
    Units = Angstrom
    Medium {
    SolverType = IEFPCM
    Solvent = Methanol
    Nonequilibrium = True
    }
    Cavity {
    Type = GePol
    Area = 1.0
    }
"""
 
# molecule
CONF_1 = psi4.geometry("""0 1
C          -0.91890         3.07720         0.54030
C          -1.74670         4.02460         1.12220
N          -1.53940         5.18890         0.40010
C          -0.61080         4.90540        -0.56980
N          -0.21640         3.64290        -0.51510
C          -0.72830         1.62840         0.90490
C          -1.26460         1.31990         2.31700
C          -1.35230         0.66500        -0.11260
C          -2.49740         1.06730        -0.82350
C          -3.13420         0.20340        -1.71730
C          -2.62480        -1.08450        -1.91310
C          -1.48570        -1.51770        -1.21710
C          -0.84220        -0.64330        -0.30220
C           0.36920        -1.16250         0.44100
C          -0.95540        -2.91350        -1.44500
H           0.36490         1.46760         0.90340
H          -2.43910         3.98610         1.95960
H          -1.98630         6.08630         0.57310
H          -0.26160         5.65570        -1.27930
H          -1.05780         0.27300         2.59250
H          -2.35950         1.45730         2.35200
H          -0.80320         1.98010         3.07170
H          -2.87940         2.08210        -0.67480
H          -4.02220         0.53510        -2.26560
H          -3.11260        -1.76940        -2.61620
H           1.15530        -1.49800        -0.26020
H           0.82180        -0.41540         1.10750
H           0.11090        -2.04140         1.06090
H           0.08340        -2.90300        -1.82440
H          -0.94010        -3.50720        -0.51200
H          -1.57380        -3.45640        -2.17790
""", name="CONF_1")
 
psi4.set_options({
    "save_jk": True,
    "pcm": True,
    "pcm__input": pcm_string
})
 
e, wfn = psi4.energy("B3LYP-D3/6-311+g(d,p)", return_wfn=True, molecule=CONF_1)
res = tdscf_excitations(wfn, states=100)

计算完毕，可以将波长（激发能）、振子强度、转子强度打印出来以便用绘图工具绘制UV与ECD图谱。 

1
2
3
4
5
6
7
8
9

import numpy as np
 
# get poles and residues to plot OPA and ECD spectra
poles = [r["EXCITATION ENERGY"] for r in res]
opa_residues = [np.linalg.norm(r["ELECTRIC DIPOLE TRANSITION MOMENT (LEN)"])**2 for r in res]
ecd_residues = [r["ROTATORY STRENGTH (LEN)"] for r in res]
print(poles)
print(opa_residues)
print(ecd_residues)

将poles(激发能量，也就是波长), opa_residues（振子强度）与ecd_residues（转子强度）三个数组打印，结果如下:

1
2
3

[0.1863171096523767, 0.1939039458434858, 0.19659912134475524, 0.2050404252248169, 0.2078159462592276, 0.21483312155800743, 0.21558247810441128, 0.21629287700878633, 0.22131073136342963, 0.2239578317414286, 0.2256655366089712, 0.2266063071494987, 0.22719511927490812, 0.22912871681783578, 0.23036097855538518, 0.23196639438645655, 0.23288762018732767, 0.23501998580463102, 0.23608956077502005, 0.23782839085338345, 0.24087960247833654, 0.24129619231009988, 0.24312235507952845, 0.24410980694671544, 0.24589867812187582, 0.24602973193576907, 0.2474711509867622, 0.24806441736255788, 0.24915399324363377, 0.2505334198822149, 0.2507468669141498, 0.2516462192853461, 0.2529149024536076, 0.2539699801159058, 0.25493797633166865, 0.25604830696224573, 0.256712318336855, 0.25910518815857936, 0.25934435981982784, 0.26037322293828696, 0.26064815699912447, 0.26262871872268084, 0.26279939036390276, 0.26427965834506684, 0.2648107191555363, 0.2651422483600908, 0.2669828436118275, 0.26821836857864967, 0.26923280469950706, 0.2700501638430005, 0.2716009041123941, 0.2726448138481998, 0.2732880554160093, 0.2736440419875829, 0.274309933278269, 0.2744231250304597, 0.2746492115721138, 0.2754770447976371, 0.27620293298686066, 0.2769771531293081, 0.27768323690960944, 0.2785993363606018, 0.2812450116133058, 0.28204969870853414, 0.2829438615368195, 0.2829959857005867, 0.28372193197946083, 0.2850708970415412, 0.28577833597122376, 0.28654236330135197, 0.2876612810257548, 0.28882846369903425, 0.2889687498366261, 0.29223770672702615, 0.2924971774299208, 0.29301078633501404, 0.2935533342680961, 0.29517983860033936, 0.2959742796563639, 0.29693311006244916, 0.2978151847572293, 0.2989882344717956, 0.29912882378353234, 0.3003095203736919, 0.3004241483808328, 0.3006209011700486, 0.3014339300694061, 0.3020126831499768, 0.3028025674303158, 0.30344091532026335, 0.3038045324249683, 0.30432506025932426, 0.30542062626395344, 0.305848061297033, 0.3065622201511003, 0.3073873287851144, 0.30803302494442214, 0.30875256479223057, 0.3091063752678208, 0.3093240808604151]
[0.008660083691178442, 0.10349305939435005, 0.048931279920771016, 0.01664861869646239, 0.023607697352347273, 0.4331760971827799, 0.19722087693559762, 0.1610869707851778, 0.5004291874968285, 0.5362120773474215, 0.7385523305745509, 1.0827967793656208, 1.0953285186803543, 1.2331583251069618, 1.7809304053128256, 0.9159936256376622, 0.15939112329899552, 1.2124717270895917, 0.08368402618706905, 0.6349603125759534, 0.16955022615511947, 0.004088688844779394, 0.11457343408462559, 0.5537920767918361, 0.008431317491542988, 0.055311354295249686, 0.09490812203280845, 0.14436463677335507, 0.12293991809471337, 0.03846757998370501, 0.04355763322776323, 0.03549218640806217, 0.018829640362198714, 0.015675566949828388, 0.06707526979885038, 0.0472314396469641, 0.001410585141789898, 0.027249991469494337, 0.01215669411084908, 0.02896998022629227, 0.023703964733904008, 0.05161528208433824, 0.06505625050737211, 0.06650963753683756, 0.3069988672801764, 0.024996907357562683, 0.013209539481416309, 0.052133675150076955, 0.049847755758164566, 0.26341816888027014, 0.007457044778791769, 0.005532852603604105, 0.0029111373742524587, 0.13346151626882374, 0.07063898473568661, 0.05371989166671687, 0.020279016086313446, 0.05534948465653796, 0.13935646199510285, 0.10236858715359076, 0.03602978482437635, 0.12279445737442111, 0.0011384301920499532, 0.1669445184936711, 0.13900267330987612, 0.003922858485403657, 0.12277084337674618, 0.0901787539114518, 0.08064482953366019, 0.0012033748240366372, 0.05832423778338989, 0.026998761335166, 0.3057183938678913, 0.15798173989573985, 0.055852891878045915, 0.17497903314275592, 0.037460505497793214, 0.09543486080372354, 0.023497561837022337, 0.04359603931558043, 0.07871790650307069, 0.050023076497197116, 0.04721549532703174, 0.001442058515559898, 0.007704308350750135, 0.032359280805689984, 0.030062327596333697, 0.1423759674141197, 0.30719405791743637, 0.08069097181483685, 0.10063493272630698, 0.09834923258949534, 0.02403182544587331, 0.18789039951589168, 0.021025667764258474, 0.03665697332342831, 0.11860141677048568, 0.013461398062255187, 0.023357460709835222, 0.008771811915421562]
[-0.0023185573599461603, 0.007772880795271462, -0.023006054965657698, 0.02823403191495033, -0.0030940906925422634, 0.033226070133156656, 0.026959457729096383, -0.007603433402276863, -0.0033942796205106507, 0.12709248350914848, 0.006601752756492565, -0.12323723569666853, 0.09016329672615413, -0.17172020363322477, 0.09223222988505309, 0.11378531083618075, -0.018570489116571665, 0.07928101788937807, -0.0362165571670183, -0.17683827015687828, -0.03409388329565155, 0.002638020701523085, 0.030585166885200642, 0.021051040289371256, 0.008822610641127486, -0.012552096713053363, 0.03546296455003074, -0.027230305756463534, -0.04345447982360856, 0.018979307688783805, -0.0205165787254387, 0.005159670113488453, -0.002044918212805822, 0.00013773234326517456, -0.028981342251422307, -0.014244631324884715, -0.001852254650470736, -0.016698127679452134, 0.010445765682826128, 5.988789893958402e-05, 0.0055470568610092585, 0.01099414941510461, -0.024311972717794596, -0.025265724796203617, 0.0037465736174804986, -0.002780122636042534, 0.008249586231930037, -0.018901303345790908, 0.013273229076012794, 0.0034774907690217015, -0.012259300424381028, -0.00038453055996658765, 0.0013796841993837599, 0.02460696753483636, 0.017809394239199774, 0.02147496433000713, 0.014033966718485457, 0.008078207634736823, 0.024635290576418802, 0.04407696637364731, -0.009629366290395465, -0.111195034016115, 0.0010305382280823652, -0.03962034048826819, 0.022109867011012085, -0.003894319272733941, -0.009897340425816765, 0.01983557425706388, 0.026573444857911882, 0.001187866266289116, 0.01631066940269573, 0.005504560000458998, -0.03100658728361818, -0.008949408647276294, 0.05188512541525664, -0.013211560802535487, -0.02509309024064607, -0.013474468990270466, 0.0058157801608803225, 0.00941595475192742, -0.001600915522050994, 0.0006496899062446836, 0.00570793877230435, -0.0015902489599031139, 0.0020201232011859624, 0.02135661646865954, -0.02734191973113685, 0.0558643389082213, -0.025157937409131477, 0.03787978186567708, -0.02756321082487943, -0.009788330784436538, 0.01812071052000217, 0.008672496271963774, -0.00023955892922729403, 0.0008040016414064197, -0.01293985186918948, 0.006059973401203619, 0.034157745289823654, -0.0025525230946705493]

poles, opa_residues与ecd_residues的信息还可以在输出文件CONF_1_tddft.out找到，部分结果如下图9所示。

用Flare QM计算化合物的ECD图谱-墨灵格的博客

图9. 激发态计算结果文件里关于激发能、振子强度与转子强度的信息，这里只是呈现了前10个跃迁。

有了激发能、振子强度与转子强度，就可以绘制UV与ECD图谱，代码如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19

import numpy as np
from rdkit import Chem
from rdkit.Chem import altair_spectrum
from altair_spectrum import plot_spectrum
from psi4.driver.p4util import spectrum
 
opa_spectrum = spectrum(poles=poles, residues=opa_residues, gamma=0.01, out_units="nm")
ecd_spectrum = spectrum(poles=poles, residues=ecd_residues, kind="ECD", gamma=0.01, out_units="nm")
 
opa_plot = plot_spectrum(opa_spectrum,
                         title="OPA (Gaussian broadening)",
                         x_title=("λ", "nm"))
 
ecd_plot = plot_spectrum(ecd_spectrum,
                         title="ECD (Gaussian broadening)",
                         x_title=("λ", "nm"),
                         y_title=("Δε", "L⋅mol⁻¹⋅cm⁻¹"))
 
opa_plot & ecd_plot

激发态计算相当快，在我的24核心机器上大约1小时不到就完成了。计算图谱如图10所示，需要注意的是：sigma值是经验值，这意味着需要调整该值以使得计算图谱与实验图谱尽可能一致。

用Flare QM计算化合物的ECD图谱-墨灵格的博客

图10. 构象1的计算ECD图谱。上：UV谱；下：ECD谱。

每个构象都有一个自己的ECD图谱，需要进一步采用玻尔兹曼加权平均法获得多个构象的平均图谱，最后将实验图谱与平均图谱进行比较以确认立体化合物构型。具体过程请参见前文[8]，这里不再叙述。

小结

总的来说，Flare集成了ECD计算所需的全部工具，可以一步到位的实现3D结构创建、XED力场结构优化、XED力场构象搜索、QM构象优化、激发态计算与UV、ECD图谱绘制。整个过程操作便利、容易实现、而且计算高效。Flare是非常实用的ECD计算辅助结构确证平台。

文献

  1. Flare V6发布——新特性、新功能与改进. http://blog.molcalx.com.cn/2022/07/06/flare-v6-released.html
  2. Flare V6预告——计算化学先睹为快. http://blog.molcalx.com.cn/2022/05/03/sneak-peek-flare-v6-comp-chem.html
  3. Flare V6预告——药化先睹为快. http://blog.molcalx.com.cn/2022/04/26/sneak-peek-flare-v6-medicinal-chemists.html
  4. 势能面扫描在Torsion Profile中的应用. http://blog.molcalx.com.cn/2020/01/30/qm-torsion-profile.html
  5. PsiCode-PSI4: Ab Initio Quantum Chemistry. https://psicode.org
  6. Roberto Di Remigio. TDDFT tutorial at PsiCon2020. https://github.com/robertodr/tddft_psicon2020
  7. Pescitelli, G.; Bruhn, T. Good Computational Practice in the Assignment of Absolute Configurations by TDDFT Calculations of ECD Spectra. Chirality 2016, 28 (6), 466–474. https://doi.org/10.1002/chir.22600.
  8. ECD图谱预测与生成Excel报告. http://blog.molcalx.com.cn/2020/01/01/ecd.html

联系我们、获取试用