CN107004064A - 自由能计算装置、方法、程序、以及记录有该程序的记录介质 - Google Patents

自由能计算装置、方法、程序、以及记录有该程序的记录介质 Download PDF

Info

Publication number
CN107004064A
CN107004064A CN201580052788.3A CN201580052788A CN107004064A CN 107004064 A CN107004064 A CN 107004064A CN 201580052788 A CN201580052788 A CN 201580052788A CN 107004064 A CN107004064 A CN 107004064A
Authority
CN
China
Prior art keywords
set body
atom
atom set
interaction energy
coordinate
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN201580052788.3A
Other languages
English (en)
Other versions
CN107004064B (zh
Inventor
松林伸幸
增田友秀
谷村隆次
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Osaka University NUC
Toray Industries Inc
Original Assignee
Osaka University NUC
Toray Industries Inc
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Osaka University NUC, Toray Industries Inc filed Critical Osaka University NUC
Publication of CN107004064A publication Critical patent/CN107004064A/zh
Application granted granted Critical
Publication of CN107004064B publication Critical patent/CN107004064B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B5/00—ICT specially adapted for modelling or simulations in systems biology, e.g. gene-regulatory networks, protein interaction networks or metabolic networks
    • G16B5/30—Dynamic-time models
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/18—Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis
    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B5/00—ICT specially adapted for modelling or simulations in systems biology, e.g. gene-regulatory networks, protein interaction networks or metabolic networks
    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B5/00—ICT specially adapted for modelling or simulations in systems biology, e.g. gene-regulatory networks, protein interaction networks or metabolic networks
    • G16B5/20—Probabilistic models
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/10—Numerical modelling

Abstract

本发明的课题在于：提供高速且高精度地计算不同分子间的自由能差的装置、方法、程序和记录有该程序的记录介质，为了实现所述课题，在本发明涉及的装置、方法、程序和记录有该程序的记录介质中，根据数学式(1)来计算ΔG。

Description

自由能计算装置、方法、程序、以及记录有该程序的记录介质
技术领域
本发明涉及一种高速且高精度地计算不同分子间的自由能差的装置、方法、程序和记录有该程序的记录介质。
背景技术
由于X射线晶体结构分析、核磁共振(NMR)等立体结构分析技术和基因工程学技术的发展，作为药物研发的靶向的蛋白的立体结构逐渐明了。随之，根据蛋白的立体结构，人们对激活或抑制该蛋白的生物体功能(生物功能)的配体的分子结构正在积极地进行开发。
作为根据蛋白的立体结构来开发配体的分子结构的特别重要的技术之一，可以列举分子动力学法或蒙特卡罗法，其是利用计算机模拟，通过计算蛋白与相对于该蛋白的不同配体间的结合自由能差来预测结合亲和性。这里，蛋白与配体间的结合自由能是与作为蛋白和配体的结合亲和性的指标的结合常数密切相关的物理量，相对于同一蛋白的不同配体间的结合自由能差是和这些不同的配体与蛋白的结合常数之比密切相关的物理量。因此，根据结合自由能差，可以求出相对于同一蛋白的不同配体间的结合亲和性之差。即，如果能够高速且高精度地计算同一蛋白与配体间的结合自由能差，则在配体的合成和结合强度实验之前，可以设计配体的分子结构，使配体的官能团或骨架结构发生变化，使蛋白与配体的自由能变小，可以期待药物开发的加速化和效率化。通常，药物开发要求一种高速且高精度的自由能计算方法，所述高速是指达到可以用数天以下的实用的计算时间进行计算的程度，而高精度是指达到能够以±1kcal/mol以下的平均误差重现通过结合强度实验获得的结合自由能或结合自由能差的程度。
作为代表性的自由能计算方法，可以列举：自由能微扰法(非专利文献1)、颗粒插入法(非专利文献2)、能量表示法(能量显示法)(非专利文献3～5)。
自由能微扰法是指如下的计算方法：逐渐改变配体与水或蛋白等配体的周边环境的分子间相互作用的相关参数，导入连接计算对象的始状态与终状态的多个假想的中间状态，分别求出相邻状态间的自由能差并进行合计，从而求出作为目标的始状态与终状态间的自由能差。通常，自由能微扰法可以高精度地重现通过结合强度实验获得的相对于同一蛋白的不同配体间的结合自由能差。但是存在着如下问题：需要超长的计算时间，不是实用的计算方法。
颗粒插入法是指如下的方法：在计算结合自由能时，在不含配体而由溶剂及蛋白等构成的始状态、与配体和蛋白结合形成复合体的终状态之间无需导入中间状态。另外，该颗粒插入法具有通过反复进行采样工序来生成多个终状态的工序，所述采样工序是指在始状态下随机选择位置坐标，并在该位置坐标上添加配体。然而，配体通常是在蛋白的特定活性部位与构成蛋白的特定氨基酸残基发生相互作用，配体与蛋白结合而形成复合体。因此，通过在随机选择的位置坐标上添加配体来生成配体与蛋白结合的复合体的概率非常低，因此现实中不可能生成多个该复合体，无法计算结合自由能。
能量表示法是指如下的方法：与自由能微扰法相比虽然计算精度相同，但可以高速进行自由能计算。然而，在计算结合自由能时，在其计算过程中，与颗粒插入法一样，具有通过反复进行采样工序来生成多个由配体与蛋白结合的复合体和溶剂构成的终状态的工序，所述采样工序：是指在不含配体而由溶剂及蛋白等构成的始状态下随机选择位置坐标、并在该位置坐标上添加配体，因此该方法无法计算结合自由能。
现有技术文献
非专利文献
非专利文献1：The Journal of Chemical Physics，22(8)，1420-1426 (1954)；
非专利文献2：The Journal of Chemical Physics, 39(11), 2808-2812 (1963)；
非专利文献3：The Journal of Chemical Physics, 113(15), 6070-6081 (2000)；
非专利文献4：The Journal of Chemical Physics, 117(8), 3605-3616 (2002)；
非专利文献5：The Journal of Chemical Physics, 119(18), 9686-9702 (2003)。
发明内容
发明所要解决的课题
在以往的自由能计算方法中，难以高速且高精度地计算结合自由能。特别是，自由能微扰法需要导入多个连结始状态和终状态的中间状态，存在着计算时间变长的问题。另外，在无需导入中间状态的颗粒插入法及能量表示法中，由于要在不含配体而由溶剂及蛋白等构成的始状态下随机选择的位置坐标上添加配体，所以无法生成多个配体与蛋白结合的复合体，存在着无法计算结合自由能的问题。因此，本发明的目的在于：提供高速且高精度地计算不同分子间的自由能差的装置、方法、程序和记录有该程序的记录介质。
用于解决课题的手段
为了解决上述课题，本发明提供一种计算装置，其具备控制部，关于反应式(1)所表示的变化，上述控制部计算变化前的原子集合体A的自由能和片段B的自由能的和、与变化后的原子集合体AB的自由能之差ΔG，
[式中，A表示由结构a构成或包含结构a而构成的原子集合体，B表示由1个以上的原子构成的片段，AB表示由原子集合体A和与该原子集合体A的结构a结合的片段B构成的原子集合体。]
其中，控制部具备下述单元：
第1原子集合体模型制作单元，制作将变化前的原子集合体A进行模型化而获得的第1原子集合体模型；
第1坐标获取单元，通过针对由第1原子集合体模型制作单元制作的第1原子集合体模型的计算机模拟，获取第1～第i的状态F1～Fi (i为2以上的整数。)的各状态下的原子集合体A的坐标；
第2坐标获取单元，根据通过第1坐标获取单元获取的原子集合体A的坐标，获取通过片段B与状态F1～Fi的各状态下的原子集合体A结合而产生的原子集合体AB的坐标；
第1相互作用能φ频率分布制作单元，根据通过第2坐标获取单元获取的原子集合体AB的坐标，计算结构a和与该结构a结合的片段B的相互作用能φ，制作表示相互作用能φ的各级频率的频率分布；
第1相互作用能φ出现概率计算单元，根据通过第1相互作用能φ频率分布制作单元制作的频率分布，计算相互作用能φ的各级的出现概率P0(φ)；
第1相互作用能ε频率分布制作单元，根据通过第2坐标获取单元获取的原子集合体AB的坐标，计算通过从原子集合体AB中去除由结构a和与该结构a结合的片段B构成的结构aB而产生的原子集合体的一部分或全部与片段B的相互作用能ε，制作通过第1相互作用能φ频率分布制作单元制作的频率分布的相互作用能φ的各级中的、表示相互作用能ε的各级频率的频率分布；
第1相互作用能ε出现概率计算单元，根据通过第1相互作用能ε频率分布制作单元制作的频率分布，计算在相互作用能φ的各级中相互作用能ε的各级的出现概率P0’(ε；φ)；
第2原子集合体模型制作单元，制作将变化后的原子集合体AB模型化而获得的第2原子集合体模型；
第3坐标获取单元，通过针对由第2原子集合体模型制作单元制作的第2原子集合体模型的计算机模拟，获取第1～第j的状态G1～Gj (j为2以上的整数。)的各状态下的原子集合体AB的坐标；
第2相互作用能φ频率分布制作单元，根据通过第3坐标获取单元获取的原子集合体AB的坐标，计算结构a和与该结构a结合的片段B的相互作用能φ，制作表示相互作用能φ的各级频率的频率分布；
第2相互作用能φ出现概率计算单元，根据通过第2相互作用能φ频率分布制作单元制作的频率分布，计算相互作用能φ的各级的出现概率P(φ)；
第2相互作用能ε频率分布制作单元，根据通过第3坐标获取单元获取的原子集合体AB的坐标，计算通过从原子集合体AB中去除由结构a和与该结构a结合的片段B构成的结构aB而产生的原子集合体的一部分或全部与片段B的相互作用能ε，制作表示相互作用能ε的各级频率的频率分布；
第2相互作用能ε出现概率计算单元，根据通过第2相互作用能ε频率分布制作单元制作的频率分布，计算相互作用能ε的各级的出现概率P’(ε)；
∫Δν(φ)P(φ)dφ计算单元，根据通过第2相互作用能φ出现概率计算单元计算的P(φ)、通过第1相互作用能ε出现概率计算单元计算的P0’(ε；φ)、以及通过第2相互作用能ε出现概率计算单元计算的P’(ε)，计算相互作用能ε所引起的自由能变化量∫Δν(φ)P(φ)dφ[式中，Δν(φ)表示相互作用能φ的各级中的相互作用能ε所引起的自由能变化量。]；以及
ΔG计算单元，根据通过第1相互作用能φ出现概率计算单元计算的P0(φ)、通过第2相互作用能φ出现概率计算单元计算的P(φ)、通过∫Δν(φ)P(φ)dφ计算单元计算的∫Δν(φ)P(φ)dφ和数学式(1)，来计算ΔG：
[数学式1]
[式中，R表示气体常数，T表示反应式(1)所表示的变化发生时的绝对温度。]。
在本发明的计算装置的一个方案(以下称作“方案1A”)中，
第2坐标获取单元进行如下事项：
制作将原子集合体C进行模型化而获得的第3原子集合体模型，所述原子集合体C由结构a和与该结构a结合的片段B构成或者包含结构a和与该结构a结合的片段B而构成；
通过针对所制作的第3原子集合体模型的计算机模拟，获取第1～第k的状态H1～Hk (k为2以上的整数。)的各状态下的原子集合体C的坐标；
选择由选自构成结构a的原子的1个或2个以上的原子构成的选择原子组，使选自状态H1～Hk的1或2个以上的状态下的原子集合体C的选择原子组的坐标相对于状态F1～Fi的各状态下的原子集合体A的选择原子组的坐标旋转和/或并进，从而在原子集合体A的选择原子组与原子集合体C的选择原子组之间制作对应的原子间距离的平方和达到最小的原子集合体C的坐标，根据所制作的原子集合体C的坐标，使选自状态H1～Hk的1个或2个以上的状态下的原子集合体C与原子集合体A重叠，根据原子集合体A的坐标和与原子集合体A重叠的原子集合体C的片段B的1个或2个以上的坐标，获取通过片段B与原子集合体A结合而产生的1个或2个以上的原子集合体AB的坐标。
在本发明的计算装置的一个方案(以下称作“方案2A”)中，
∫Δν(φ)P(φ)dφ计算单元根据通过第2相互作用能φ出现概率计算单元计算的P(φ)、通过第1相互作用能ε出现概率计算单元计算的P0’(ε；φ)、以及通过第2相互作用能ε出现概率计算单元计算的P’(ε)，利用能量表示法计算相互作用能ε所引起的自由能变化量∫Δν(φ)P(φ)dφ。
在本发明的计算装置的一个方案(以下称作“方案3A”)中，
片段B由包含作为虚拟原子的点电荷的原子构成，第2坐标获取单元在构成原子集合体A的结构a的原子的电荷参数上添加片段B的点电荷。
在本发明的计算装置中，可以将方案1A～方案3A中的2个以上的方案组合。
另外，本发明还提供一种计算方法，关于反应式(1)所表示的变化，上述计算方法利用计算机计算变化前的原子集合体A的自由能和片段B的自由能的和、与变化后的原子集合体AB的自由能之差ΔG，
[式中，A表示由结构a构成或包含结构a而构成的原子集合体，B表示由1个以上的原子构成的片段，AB表示由原子集合体A和与该原子集合体A的结构a结合的片段B构成的原子集合体。]，
其中，计算机的控制部实行如下步骤：
第1原子集合体模型制作步骤，制作将变化前的原子集合体A进行模型化而获得的第1原子集合体模型；
第1坐标获取步骤，通过针对第1原子集合体模型制作步骤制作的第1原子集合体模型的计算机模拟，获取第1～第i的状态F1～Fi (i为2以上的整数。)的各状态下的原子集合体A的坐标；
第2坐标获取步骤，根据通过第1坐标获取步骤获取的原子集合体A的坐标，获取通过片段B与状态F1～Fi的各状态下的原子集合体A结合而产生的原子集合体AB的坐标；
第1相互作用能φ频率分布制作步骤，根据通过第2坐标获取步骤获取的原子集合体AB的坐标，计算结构a和与该结构a结合的片段B的相互作用能φ，制作表示相互作用能φ的各级频率的频率分布；
第1相互作用能φ出现概率计算步骤，根据通过第1相互作用能φ频率分布制作步骤制作的频率分布，计算相互作用能φ的各级的出现概率P0(φ)；
第1相互作用能ε频率分布制作步骤，根据通过第2坐标获取步骤获取的原子集合体AB的坐标，计算通过从原子集合体AB中去除由结构a和与该结构a结合的片段B构成的结构aB而产生的原子集合体的一部分或全部与片段B的相互作用能ε，制作通过第1相互作用能φ频率分布制作步骤制作的频率分布的相互作用能φ的各级中的、表示相互作用能ε的各级频率的频率分布；
第1相互作用能ε出现概率计算步骤，根据通过第1相互作用能ε频率分布制作步骤制作的频率分布，计算相互作用能φ的各级中的相互作用能ε的各级的出现概率P0’(ε；φ)；
第2原子集合体模型制作步骤，制作将变化后的原子集合体AB模型化而获得的第2原子集合体模型；
第3坐标获取步骤，通过针对第2原子集合体模型制作步骤中制作的第2原子集合体模型的计算机模拟，获取第1～第j的状态G1～Gj (j为2以上的整数。)的各状态下的原子集合体AB的坐标；
第2相互作用能φ频率分布制作步骤，根据通过第3坐标获取步骤获取的原子集合体AB的坐标，计算结构a和与该结构a结合的片段B的相互作用能φ，制作表示相互作用能φ的各级频率的频率分布；
第2相互作用能φ出现概率计算步骤，根据通过第2相互作用能φ频率分布制作步骤制作的频率分布，计算相互作用能φ的各级的出现概率P(φ)；
第2相互作用能ε频率分布制作步骤，根据通过第3坐标获取步骤获取的原子集合体AB的坐标，计算通过从原子集合体AB中去除由结构a和与该结构a结合的片段B构成的结构aB而产生的原子集合体的一部分或全部与片段B的相互作用能ε，制作表示相互作用能ε的各级频率的频率分布；
第2相互作用能ε出现概率计算步骤，根据通过第2相互作用能ε频率分布制作步骤制作的频率分布，计算相互作用能ε的各级的出现概率P’(ε)；
∫Δν(φ)P(φ)dφ计算步骤，根据通过第2相互作用能φ出现概率计算步骤计算的P(φ)、通过第1相互作用能ε出现概率计算步骤计算的P0’(ε；φ)、以及通过第2相互作用能ε出现概率计算步骤计算的P’(ε)，计算相互作用能ε所引起的自由能变化量∫Δν(φ)P(φ)dφ [式中，Δν(φ)表示相互作用能φ的各级中的、相互作用能ε所引起的自由能变化量。]；
ΔG计算步骤，根据通过第1相互作用能φ出现概率计算步骤计算的P0(φ)、通过第2相互作用能φ出现概率计算步骤计算的P(φ)、通过∫Δν(φ)P(φ)dφ计算步骤计算的∫Δν(φ)P(φ)dφ和数学式(1)，来计算ΔG：
[数学式2]
[式中，R表示气体常数，T表示反应式(1)所表示的变化发生时的绝对温度。]。
在本发明的计算方法的一个方案(以下称作“方案1B”)中，
在第2坐标获取步骤中，计算机的控制部实行如下事项：
制作将原子集合体C进行模型化而获得的第3原子集合体模型，上述原子集合体C由结构a和与该结构a结合的片段B构成或者包含结构a和与该结构a结合的片段B而构成；
通过针对所制作的第3原子集合体模型的计算机模拟，获取第1～第k的状态H1～Hk (k为2以上的整数。)的各状态下的原子集合体C的坐标；
选择由选自构成结构a的原子的1个或2个以上的原子构成的选择原子组，使选自状态H1～Hk的1种或2种以上的状态下的原子集合体C的选择原子组的坐标相对于状态F1～Fi的各状态下的原子集合体A的选择原子组的坐标旋转和/或并进，从而在原子集合体A的选择原子组与原子集合体C的选择原子组之间制作对应的原子间距离的平方和达到最小的原子集合体C的坐标，根据所制作的原子集合体C的坐标，使选自状态H1～Hk的1种或2种以上的状态下的原子集合体C与原子集合体A重叠，根据原子集合体A的坐标和与原子集合体A重叠的原子集合体C的片段B的1个或2个以上的坐标，获取通过片段B与原子集合体A结合而产生的1个或2个以上的原子集合体AB的坐标。
在本发明的计算方法的一个方案(以下称作“方案2B”)中，
在∫Δν(φ)P(φ)dφ计算步骤中，计算机的控制部根据通过第2相互作用能φ出现概率计算步骤计算的P(φ)、通过第1相互作用能ε出现概率计算步骤计算的P0’(ε；φ)、以及通过第2相互作用能ε出现概率计算步骤计算的P’(ε)，利用能量表示法计算相互作用能ε所引起的自由能变化量∫Δν(φ)P(φ)dφ。
在本发明的计算方法的一个方案(以下称作“方案3B”)中，
片段B由包含作为虚拟原子的点电荷的原子构成，在第2坐标获取步骤中，计算机的控制部将片段B的点电荷添加在构成原子集合体A的结构a的原子的电荷参数中。
在本发明的计算方法中，可以将方案1B～方案3B中的2个以上的方案组合。
而且，本发明还提供一种程序，关于反应式(1)所表示的变化，上述程序计算变化前的原子集合体A的自由能和片段B的自由能的和、与变化后的原子集合体AB的自由能之差ΔG，
[式中，A表示由结构a构成或者包含结构a而构成的原子集合体，B表示由1个以上的原子构成的片段，AB表示由原子集合体A和与该原子集合体A的结构a结合的片段B构成的原子集合体。]，
上述程序使计算机的控制部起到如下单元的作用：
第1原子集合体模型制作单元，制作将变化前的原子集合体A进行模型化而获得的第1原子集合体模型；
第1坐标获取单元，通过针对由第1原子集合体模型制作单元制作的第1原子集合体模型的计算机模拟，获取第1～第i的状态F1～Fi (i为2以上的整数。)的各状态下的原子集合体A的坐标；
第2坐标获取单元，根据通过第1坐标获取单元获取的原子集合体A的坐标，获取通过片段B与状态F1～Fi的各状态下的原子集合体A结合而产生的原子集合体AB的坐标；
第1相互作用能φ频率分布制作单元，根据通过第2坐标获取单元获取的原子集合体AB的坐标，计算结构a和与该结构a结合的片段B的相互作用能φ，制作表示相互作用能φ的各级频率的频率分布；
第1相互作用能φ出现概率计算单元，根据通过第1相互作用能φ频率分布制作单元制作的频率分布，计算相互作用能φ的各级的出现概率P0(φ)；
第1相互作用能ε频率分布制作单元，根据通过第2坐标获取单元获取的原子集合体AB的坐标，计算通过从原子集合体AB中去除由结构a和与该结构a结合的片段B构成的结构aB而产生的原子集合体的一部分或全部与片段B的相互作用能ε，制作通过第1相互作用能φ频率分布制作单元制作的频率分布的相互作用能φ的各级中的、表示相互作用能ε的各级频率的频率分布；
第1相互作用能ε出现概率计算单元，根据通过第1相互作用能ε频率分布制作单元制作的频率分布，计算在相互作用能φ的各级中相互作用能ε的各级的出现概率P0’(ε；φ)；
第2原子集合体模型制作单元，制作将变化后的原子集合体AB进行模型化而获得的第2原子集合体模型；
第3坐标获取单元，通过针对由第2原子集合体模型制作单元制作的第2原子集合体模型的计算机模拟，获取第1～第j的状态G1～Gj (j为2以上的整数。)的各状态下的原子集合体AB的坐标；
第2相互作用能φ频率分布制作单元，根据通过第3坐标获取单元获取的原子集合体AB的坐标，计算结构a和与该结构a结合的片段B的相互作用能φ，制作表示相互作用能φ的各级频率的频率分布；
第2相互作用能φ出现概率计算单元，根据通过第2相互作用能φ频率分布制作单元制作的频率分布，计算相互作用能φ的各级的出现概率P(φ)；
第2相互作用能ε频率分布制作单元，根据通过第3坐标获取单元获取的原子集合体AB的坐标，计算通过从原子集合体AB中去除由结构a和与该结构a结合的片段B构成的结构aB而产生的原子集合体的一部分或全部与片段B的相互作用能ε，制作表示相互作用能ε的各级频率的频率分布；
第2相互作用能ε出现概率计算单元，根据通过第2相互作用能ε频率分布制作单元制作的频率分布，计算相互作用能ε的各级的出现概率P’(ε)；
∫Δν(φ)P(φ)dφ计算单元，根据通过第2相互作用能φ出现概率计算单元计算的P(φ)、通过第1相互作用能ε出现概率计算单元计算的P0’(ε；φ)、以及通过第2相互作用能ε出现概率计算单元计算的P’(ε)，计算相互作用能ε所引起的自由能变化量∫Δν(φ)P(φ)dφ [式中，Δν(φ)表示相互作用能φ的各级中的、相互作用能ε所引起的自由能变化量。]；以及
ΔG计算单元，根据通过第1相互作用能φ出现概率计算单元计算的P0(φ)、通过第2相互作用能φ出现概率计算单元计算的P(φ)、通过∫Δν(φ)P(φ)dφ计算单元计算的∫Δν(φ)P(φ)dφ和数学式(1)，来计算ΔG：
[数学式3]
[式中，R表示气体常数，T表示反应式(1)所表示的变化发生时的绝对温度。]。
在本发明的程序的一个方案(以下称作“方案1C”)中，
第2坐标获取单元实行如下事项：
制作将由结构a和与该结构a结合的片段B构成或者包含结构a和与该结构a结合的片段B而构成的原子集合体C进行模型化而获得的第3原子集合体模型；
通过针对所制作的第3原子集合体模型的计算机模拟，获取第1～第k的状态H1～Hk (k为2以上的整数。)的各状态下的原子集合体C的坐标；
选择由选自构成结构a的原子的1个或2个以上的原子构成的选择原子组，使选自状态H1～Hk的1种或2种以上的状态下的原子集合体C的选择原子组的坐标相对于状态F1～Fi的各状态下的原子集合体A的选择原子组的坐标旋转和/或并进，从而在原子集合体A的选择原子组与原子集合体C的选择原子组之间制作对应的原子间距离的平方和达到最小的原子集合体C的坐标，根据所制作的原子集合体C的坐标，使选自状态H1～Hk的1种或2种以上的状态下的原子集合体C与原子集合体A重叠，根据原子集合体A的坐标和与原子集合体A重叠的原子集合体C的片段B的1个或2个以上的坐标，获取通过片段B与原子集合体A结合而产生的1个或2个以上的原子集合体AB的坐标。
在本发明的程序的一个方案(以下称作“方案2C”)中，
∫Δν(φ)P(φ)dφ计算单元根据通过第2相互作用能φ出现概率计算单元计算的P(φ)、通过第1相互作用能ε出现概率计算单元计算的P0’(ε；φ)、以及通过第2相互作用能ε出现概率计算单元计算的P’(ε)，利用能量表示法计算相互作用能ε所引起的自由能变化量∫Δν(φ)P(φ)dφ。
在本发明的程序的一个方案(以下称作“方案3C”)中，
片段B由包含作为虚拟原子的点电荷的原子构成，第2坐标获取单元将片段B的点电荷添加在构成原子集合体A的结构a的原子的电荷参数中。
在本发明的程序中，可以将方案1C～方案3C中的2个以上的方案组合。
而且，本发明还提供一种记录有本发明的程序的计算机可读取的记录介质。
发明效果
根据本发明，通过根据原子集合体A的坐标获取因片段B与原子集合体A结合而产生的原子集合体AB的坐标，可以有效获取原子集合体AB的坐标。由此，与颗粒插入法及能量表示法相比，可以使从通过计算机模拟而产生的系综(Ensemble)中的抽样效率飞跃性地提高。而且，通过提高抽样效率，在统计学上能够以准确的精度计算满足古典(经典)统计力学理论的数学式(1)，可以确保高的计算精度。并且，由于不需要连接始状态和终状态的中间状态，所以与自由能微扰法相比可以大幅缩短计算时间。因此，根据本发明，可以高速且高精度地计算不同分子间的结合自由能差，可以实现新型药品研发的加速化和效率化。
附图说明
[图1] 图1是显示本发明的一实施方式所涉及的计算装置的构成的功能方框图。
[图2] 图2是显示同一实施方式所涉及的计算装置的处理次序(处理步骤)的流程图。
[图3] 图3是显示步骤S120中的处理次序的一实施方式的流程图。
[图4] 图4是显示步骤S130中的处理次序的一实施方式的流程图。
[图5] 图5是显示步骤S130中的处理次序的另一实施方式的流程图。
[图6] 图6是显示步骤S220中的处理次序的一实施方式的流程图。
具体实施方式
通常，在古典(经典)统计力学理论中，关于反应式(1)所表示的变化，变化前的原子集合体A的自由能和片段B的自由能的和、与变化后的原子集合体AB的自由能之差用数学式(2)表示。
这里，变化前的原子集合体A和片段B处于热平衡状态，变化后的原子集合体AB也处于热平衡状态。另外，在以下的数学式(2)～(5)中，积分符号的积分(∫)是集合表示与各积分变数相关的积分符号。
另一方面，变化后的原子集合体AB中的结构a和与该结构a结合的片段B的相互作用能φ的出现概率(即能量分布)即P(φ)用数学式(3)表示，同时因片段B与变化前的原子集合体A结合而产生的原子集合体AB中的结构a和与该结构a结合的片段B的相互作用能φ的出现概率(即能量分布)即P0(φ)用数学式(4)表示，它们的比可以整理成数学式(5)。数学式(1)可以通过将数学式(5)变形而导出。这里，数学式(3)和(4)中的：
[数学式4]
是狄拉克(Dirac)δ函数，其是当φ与下式相等时为1、其它情形为零的函数。
[数学式5]
[数学式6]
ψ：构成片段B的原子的坐标的集合；
K：构成结构a的原子的坐标的集合；
X：构成通过从变化后的原子集合体AB中去除由结构a和与该结构a结合的片段B构成的结构aB而产生的原子集合体(即变化后的原子集合体AB中的结构aB的周边环境)的所有原子的坐标的集合、或者构成通过从变化前的原子集合体A中去除结构a而产生的原子集合体(即变化前的原子集合体A中的结构a的周边环境)的所有原子的坐标的集合；
x1：通过从变化后的原子集合体AB中去除由结构a和与该结构a结合的片段B构成的结构aB而产生的原子集合体(即变化后的原子集合体AB中的结构aB的周边环境)、或者通过从变化前的原子集合体A中去除结构a而产生的原子集合体(即变化前的原子集合体A中的结构a的周边环境)中的构成原子集合体的第i个分子的坐标的集合(例如，原子集合体AB包含多个水分子，对其各水分子赋予编号时，构成第i个水分子的原子的坐标的集合)；
Ξ：构成由结构a和与该结构a结合的片段B构成的结构aB的原子间的相互作用能；
Ψ：构成片段B的原子间的相互作用能；
W：构成结构a的原子间的相互作用能；
U：构成通过从变化后的原子集合体AB中去除由结构a和与该结构a结合的片段B构成的结构aB而产生的原子集合体(即变化后的原子集合体AB中的结构aB的周边环境)的原子间、或者构成通过从变化前的原子集合体A中去除结构a而产生的原子集合体(即变化前的原子集合体A中的结构a的周边环境)的原子间的相互作用能；
v：变化后的原子集合体AB中的、由结构a和与该结构a结合的片段B构成的结构aB与该结构aB的周边环境(即通过从原子集合体AB中去除该结构aB而产生的原子集合体)的相互作用能；
u：变化前的原子集合体A中的、结构a与该结构a的周边环境(即通过从变化前的原子集合体A中去除该结构a而产生的原子集合体)的相互作用能；
R：气体常数；
T：反应式(1)所表示的变化发生时的绝对温度；
这里，Ξ、Ψ、W、U、v、u中的“相互作用能”是指与键长、键角、扭转角等所关联的原子间的键合有关的相互作用所引起的能量、和范德华相互作用或静电相互作用等非键合相互作用所引起的能量。
[数学式7]
数学式(1)的第1项表示变化后的原子集合体AB中的结构a和与该结构a结合的片段B的相互作用能φ的平均值。
数学式(1)的第1项和第2项表示反应式(1)所表示的变化中的、因在片段B和原子集合体A的结构a之间产生键合而引起的自由能变化量。
数学式(1)的第3项的Δν(φ)是指当结构a和与该结构a结合的片段B的相互作用能为φ时的、反应式(1)所表示的变化中的自由能变化量，表示片段B和从原子集合体AB中去除由结构a和与该结构a结合的片段B构成的结构aB而获得的原子集合体的一部分或全部的相互作用所引起的自由能变化量。
关于数学式(1)的第2项，变化后的原子集合体AB中的结构a和与该结构a结合的片段B的相互作用能φ的出现概率即P(φ)和通过片段B与变化前的原子集合体A结合而产生的原子集合体AB中的结构a和与该结构a结合的片段B的相互作用能φ的出现概率即P0(φ)要求P(φ)与P0(φ)之比P(φ)/P0(φ)不发散。即，在P(φ)为0＜P(φ)≦1的范围具有出现概率的φ中，P0(φ)也必需具有0＜P0(φ)≦1的出现概率。
数学式(1)的第3项的∫Δν(φ)P(φ)dφ表示反应式(1)所表示的变化中的、相互作用能ε所引起的自由能变化量，可以利用以下的数学式(6)，采用以往的能量表示法进行计算。另外，数学式(6)左边的符号与以往的能量表示法中的符号一样，表示分布函数。这里，在本发明中，当片段B包含后述的点电荷等虚拟原子时，从原子集合体AB中去除片段B是指，例如当原子集合体AB为苯甲醚(C6H5OCH3)、而片段B为甲氧基(-OCH3)和点电荷(持有该甲氧基所结合的苯环的碳原子所具有的电荷参数的虚拟原子)时，从苯甲醚中去除甲氧基和该点电荷所持有的电荷参数。即，原子集合体A是用C6H5表示的原子集合体，其是在原子集合体AB中结合有甲氧基的苯环上的碳原子所具有的电荷参数为零的原子集合体。另外，片段B与原子集合体A结合是指，例如，当片段B为甲氧基(-OCH3)和点电荷(持有该甲氧基所结合的苯环的碳原子所具有的电荷参数的虚拟原子)、而原子集合体A为C6H5所表示的原子集合体、且片段B的甲氧基所结合的苯环上的碳原子所具有的部分电荷为零时，片段B的甲氧基的氧原子与苯环上的碳原子形成共价键，而且将片段B的点电荷所持有的部分电荷(电荷参数)添加在甲氧基的氧原子所结合的苯环的碳原子的电荷参数中。需要说明的是，在本发明中，有时将原子所具有的电荷参数称作部分电荷。另外，例如有时还将从构成原子集合体AB的原子的电荷参数中去除片段B的点电荷所持有的电荷参数称作从原子集合体AB中去除片段B的点电荷。同样，将片段B的点电荷所持有的电荷参数添加在构成原子集合体AB的原子的电荷参数中有时还被称作将片段B的点电荷添加在原子集合体AB中或者添加在构成原子集合体AB的原子的电荷参数中。
数学式(1)的第3项的∫Δν(φ)P(φ)dφ可以根据：变化后的原子集合体AB中的相互作用能φ的各级的出现概率即P(φ)、在通过片段B与变化前的原子集合体A结合而产生的原子集合体AB中的相互作用能φ的各级中相互作用能ε的各级的出现概率即P0’(ε；φ)、以及变化后的原子集合体AB中的相互作用能ε的各级的出现概率即P’(ε)来计算，而不用计算Δν(φ)。此外，∫Δν(φ)P(φ)dφ例如还可如下计算：计算相互作用能φ的各级中的、相互作用能ε所引起的自由能变化量Δν(φ)，利用所计算的Δν(φ)和变化后的原子集合体AB中的相互作用能φ的各级的出现概率即P(φ)来进行计算。
[数学式8]
在本发明中，如下所详述，通过在构成原子集合体A的原子上添加构成片段B的原子，可以使抽样效率飞跃性地提高。因此，对计算Δν(φ)的方法没有特别限定，例如可以采用以往的颗粒插入法或能量表示法。需要说明的是，通过本发明计算的数学式(1)的自由能差ΔG可适用于溶剂化自由能差或结合自由能差等任意的自由能差的计算。
本发明中计算的ΔG是指与反应式(1)所表示的变化有关的自由能变化量即变化前的原子集合体A的自由能和片段B的自由能的和G1、与变化后的原子集合体AB的自由能G2之差(G2-G1)。需要说明的是，在变化前，原子集合体A与片段B没有发生相互作用。
[式中，A表示由结构a构成或者包含结构a而构成的原子集合体，B表示由1个以上的原子构成的片段，AB表示由原子集合体A和与该原子集合体A的结构a结合的片段B构成的原子集合体。]
原子集合体A和原子集合体AB是由1个或2个以上的原子构成的原子的集合体。在原子集合体A和原子集合体AB中，可以包含未与其它原子结合的1个或2个以上的原子。这里，原子也包括离子。另外，还包括虚拟原子(虚构的原子)。作为构成原子集合体A和原子集合体AB的原子间键，例如可以列举：共价键、配位键、氢键、静电相互作用、疏水性相互作用等。原子集合体A和原子集合体AB可以由1种原子间键(例如共价键)构成，也可以由2种以上的原子间键的组合(例如共价键与1种或2种以上的其它原子间键的组合)构成。
原子集合体A和原子集合体AB可以由1种或2种以上的分子构成。
作为分子，例如可以列举：配体、该配体所结合的蛋白、水分子等溶剂分子等。当原子集合体A由1种分子构成时，原子集合体AB也由1种分子构成，当原子集合体A由2种以上的分子构成时，原子集合体AB也由2种以上的分子构成。作为原子集合体A和原子集合体AB由2种以上的分子构成的情形，例如可以列举：原子集合体A和原子集合体AB除包含配体以外、还包含该配体所结合的蛋白和/或水分子等溶剂分子的情形等。当原子集合体A和原子集合体AB除包含配体以外、还包含该配体所结合的蛋白时，配体与蛋白可以形成复合体。作为该复合体中的配体与蛋白间的结合方式，例如可以列举：配位键、氢键、静电相互作用、共价键、疏水性相互作用等。当原子集合体A和原子集合体AB包含配体或配体与蛋白的复合体时，作为配体或复合体与水分子等溶剂分子间的结合方式，例如可以列举：配位键、氢键、静电相互作用、共价键、疏水性相互作用等。
原子集合体A由结构a构成、或者包含结构a而构成。当原子集合体A由1个分子构成时，结构a可以是该分子的结构整体(此时，原子集合体A由结构a构成)，也可以是该分子的一部分结构(此时，原子集合体A包含结构a而构成)。当原子集合体A由2种以上的分子构成时，结构a可以是该2种以上的分子的结构整体(此时，原子集合体A由结构a构成)，也可以是该2种以上的分子的一部分结构(此时，原子集合体A包含结构a而构成)。当原子集合体A包含蛋白与配体的复合体时，作为该复合体的一部分结构，例如可以列举：配体的结构整体或一部分、蛋白的结构整体或一部分、由配体整体或一部分与一部分蛋白构成的部分、由受体整体或一部分与一部分配体构成的部分等。
原子集合体A相当于从原子集合体AB中去除了片段B而获得的原子集合体。例如，当原子集合体AB为苯甲醚(C6H5OCH3)、而片段B为甲氧基(-OCH3)时，原子集合体A为C6H5。另外，例如当原子集合体AB为苯甲醚(C6H5OCH3)、而片段B为甲氧基(-OCH3)和点电荷(持有该甲氧基所结合的苯环的碳原子所具有的电荷参数的虚拟原子)时，原子集合体A可通过从苯甲醚中去除甲氧基和该点电荷而获得。即，原子集合体A是用C6H5表示的原子集合体，其是在原子集合体AB中结合有甲氧基的苯环上的碳原子所具有的电荷参数为零的原子集合体。
原子集合体AB由原子集合体A和与该原子集合体A的结构a结合的片段B构成。原子集合体AB相当于片段B与原子集合体A结合而获得的原子集合体。例如，当原子集合体A为从苯甲醚中去除甲氧基和点电荷(该甲氧基所结合的苯环的碳原子所具有的部分电荷)而获得的结构、而片段B由甲氧基和该点电荷构成时，原子集合体AB为苯甲醚。
构成结构a的原子数只要是1个以上即可，没有特别限定，但优选为3个以上、进一步优选为5个以上。当构成结构a的原子数为3个以上时，可有效实施片段B与原子集合体A的结构a的结合，因此可提高抽样效率。另外，对构成结构a的原子数的上限值没有特别限定。
对构成结构a的原子的种类没有特别限定。作为构成结构a的原子，例如可以列举：碳原子、氢原子、氮原子、氧原子等。结构a可以由1种原子构成，也可以由2种以上的原子构成。
对构成结构a的原子间的结合方式没有特别限定。作为构成结构a的原子间的结合方式，例如可以列举：共价键、配位键、氢键、静电相互作用、疏水性相互作用等。结构a可以由1种原子间键(例如共价键)构成，也可以由2种以上的键的组合(例如共价键与1种或2种以上的其它键的组合)构成。
在构成结构a的原子中，可以包含1个或2个以上的虚拟原子(实际不存在的原子)。结构a可以仅由实际存在原子构成，也可以仅由虚拟原子构成，还可以由实际存在原子与虚拟原子的组合构成。作为虚拟原子，例如可以列举原子价壳在形式上不具有8个电子的原子。需要说明的是，已知符合属于元素周期表的1～2族和13～18族的典型元素的原子通常形成其原子价壳在形式上具有8个电子的电子构型而满足八隅体规则，在化学上稳定存在。作为其它虚拟原子，例如可以列举：在其与其它原子之间不具有与键长、键角、扭转角等所关联的键有关的相互作用或范德华相互作用(即范德华势能为零)、但具有与其它原子间的静电相互作用(即持有电荷参数、具有库仑势能)的点电荷；不具有与其它原子间的静电相互作用(即库仑势能为零)、但在其与其它原子间具有与键长、键角、扭转角等所引起的键有关的相互作用和范德华相互作用的原子；以离子状态存在的原子；在其它原子间不具有与键长、键角、扭转角等所引起的键有关的相互作用、范德华相互作用和静电相互作用中的任一种作用的虚拟原子等。需要说明的是，以下将由范德华相互作用所引起的能量称为范德华势能、将由静电相互作用所引起的能量称为库仑势能。
构成片段B的原子数只要是1个以上即可，没有特别限定。对构成片段B的原子数的上限没有特别限定，通常为25个、优选为10个。
对构成片段B的原子的种类没有特别限定。作为构成片段B的原子，例如可以列举：碳原子、氢原子、氮原子、氧原子等。片段B可以由1种原子构成，也可以由2种以上的原子构成。
对构成片段B的原子间的结合方式没有特别限定。作为构成片段B的原子间的结合方式，例如可以列举：共价键、配位键、氢键、静电相互作用、疏水性相互作用等。片段B可以由1种原子间键(例如共价键)构成，也可以由2种以上的键的组合(例如共价键与1种或2种以上的其它键的组合)构成。
在构成片段B的原子中，可以包含1个或2个以上的虚拟原子(实际不存在的原子)。片段B可以仅由实际存在原子构成，也可以仅由虚拟原子构成，还可以由实际存在原子与虚拟原子的组合构成。作为虚拟原子，可以列举与结构a相同的具体例子。当片段B包含点电荷作为虚拟原子、且由实际存在原子与虚拟原子的组合构成时，在该原子间不具有共价键，也不具有与键长、键角、扭转角等所引起的键有关的相互作用。另外，在计算机模拟上虽然优选计算片段B的虚拟原子(点电荷)与实际存在原子间的静电相互作用，但在计算自由能的过程中将其统一时，可以计算也可以不计算。当假定在片段B的虚拟原子(点电荷)与实际存在原子之间具有共价键、且该原子间具有3个以上的共价键时，优选考虑分子动力学模拟等计算机模拟中的所谓“1-4相互作用”，计算虚拟原子(点电荷)与实际存在原子间的静电相互作用。例如，当结构aB为苯甲醚(C6H5OCH3)、而片段B为甲氧基(-OCH3)和点电荷(具有该甲氧基所结合的苯环的碳原子所具有的电荷参数的虚拟原子)时，优选计算该点电荷与甲氧基的氢原子各自间的静电相互作用。
以下，根据附图，对本发明的一实施方式所涉及的计算装置进行说明。
如图1所示，计算装置1具备输入部11、控制部12、存储部13和输出部14，它们通过系统总线相互连接。计算装置1可以通过使用通用的计算机作为基本硬件而得以实现。
输入部11例如由操作者所操作的键盘、鼠标等定点设备(指点器)构成，输入来自操作者的指示(例如程序的执行指示、处理结果的显示指示)、各处理所必需的数据的输入等各种操作信号。所输入的数据通过控制部12存储在存储部13中。
控制部12例如由CPU、RAM、ROM等构成，根据存储在存储部13中的各种数据、各种程序等计算ΔG。此时，控制部12通过执行用于计算ΔG的计算程序，起到第1原子集合体模型制作单元12A、第1坐标获取单元12B、第2坐标获取单元12C、第1相互作用能φ频率分布制作单元12D、第1相互作用能φ出现概率计算单元12E、第1相互作用能ε频率分布制作单元12F、第1相互作用能ε出现概率计算单元12G、第2原子集合体模型制作单元12H、第3坐标获取单元12I、第2相互作用能φ频率分布制作单元12J、第2相互作用能φ出现概率计算单元12K、第2相互作用能ε频率分布制作单元12L、第2相互作用能ε出现概率计算单元12M、∫Δν(φ)P(φ)dφ计算单元12N、以及ΔG计算单元12O的作用。需要说明的是，控制部12的CPU可以是如下的构成：具有多个运算芯(core)，各运算芯通过执行用于计算ΔG的计算程序而起到各种单元的作用。根据这样的构成，可以并行处理多个计算，由此可以缩短ΔG的计算所需的时间。
存储部13例如由RAM、硬盘等存储器构成，存储各种数据、各种程序等。作为存储在存储部13中的数据和程序，例如可以列举：变化前的原子集合体A的相关数据13A、第1模拟条件的相关数据13B、第1模拟程序13C、变化后的原子集合体AB的相关数据13D、第2模拟条件的相关数据13E、第2模拟程序13F等。需要说明的是，在本实施方式中，虽然第1模拟程序13C和第2模拟程序13F可以作为分开的程序来记载，但第1模拟程序13C和第2模拟程序13F可以是1个程序。
输出部14例如由显示器等构成，输出通过各种单元获取或计算的结果(例如通过ΔG计算单元12O计算的ΔG等)。
第1原子集合体模型制作单元12A制作将变化前的原子集合体A进行模型化而获得的第1原子集合体模型。通过第1原子集合体模型制作单元12A制作的第1原子集合体模型的相关数据(例如构成原子集合体A的各原子的坐标、种类、质量、部分电荷、原子间的结合信息等)通过控制部12存储在存储部13中。
第1原子集合体模型制作单元12A在制作第1原子集合体模型时，使用存储在存储部13中的变化前的原子集合体A的相关数据13A。变化前的原子集合体A的相关数据13A例如以控制部12可读取的文件形式存储在存储部13中。变化前的原子集合体A的相关数据13A只要能够制作将变化前的原子集合体A进行模型化而获得的第1原子集合体模型即可，没有特别限定。作为变化前的原子集合体A的相关数据13A，例如可以列举：构成原子集合体A的各原子的坐标、种类、质量、部分电荷、原子间的结合信息等。
第1原子集合体模型只要可以进行计算机模拟即可，没有特别限定。作为原子集合体模型，可以列举：全原子模型、珠-簧模型、联合原子模型(United atom model)等。需要说明的是，珠-簧模型是指将构成原子集合体A的单体单元视为1个珠粒(段)、并通过假想的弹簧使之结合的模型，而联合原子模型是指将氢原子包含在重原子(例如碳原子)中而作为1个原子(质点)进行处理的模型。
第1坐标获取单元12B通过以针对通过第1原子集合体模型制作单元12A制作的第1原子集合体模型的计算机模拟的结果的形式输出的抽点打印(snapshot)，获取第1～第i的状态F1～Fi (i为2以上的整数。)的各状态下的原子集合体A的坐标。这里，抽点打印包含构成原子集合体A的所有原子的坐标。即，在状态F1～Fi的各状态下的抽点打印中，分别包含构成状态F1～Fi的各状态下的原子集合体A的所有原子的坐标。以下，有时将状态F1、F2、···、Fi下的原子集合体A的坐标分别称作坐标RA(F1)、RA(F2)、···、RA(Fi)。通过第1坐标获取单元12B获取的构成原子集合体A的所有原子的坐标通过控制部12存储在存储部13中。
第1坐标获取单元12B所实行的计算机模拟只要基于统计力学理论即可，没有特别限定。作为代表性的计算机模拟，例如可以列举：分子动力学法、蒙特卡罗法、将分子动力学法或蒙特卡罗法与第一原理计算组合而获得的方法等。
第1坐标获取单元12B根据存储在存储部13中的第1原子集合体模型的相关数据、第1模拟条件的相关数据13B、第1模拟程序13C等实行计算机模拟。第1模拟条件的相关数据13B以控制部12可读取的文件形式存储在存储部13中。
模拟条件只要能够实行计算机模拟即可，没有特别限定。作为模拟条件，例如可以列举：温度条件、压力条件、势能参数的种类和值、产生的系综的种类、边界条件、抽点打印数目等的输出条件等。
作为势能参数，例如可以列举：与键长、键角、扭转角等所关联的原子间的键合有关的参数、以及与在原子间起作用的范德华相互作用和静电相互作用等有关的参数。
作为势能参数，例如可以使用：Amber、GAFF、CHARMm (注册商标)、CGenFF(CHARMm36)、DISCOVER、GROMOS、DREIDING、OPLS等公知的势能参数。
例如，在使用Amber作为势能参数时，根据针对构成原子集合体A的各原子分派的原子类型和各原子间的原子类型的配对，确定与键长、键角、范德华等有关的参数的值。
势能参数分为与原子间的结合有关的势能参数和与是否结合无关的非结合势能参数。在考虑到计算精度时，优选通过以下所述的方法确定势能参数。
与原子间的结合有关的势能参数优选通过基于第一原理计算的量子化学计算来确定，作为量子化学计算，优选采用Hartree-Fock法(以下称作“HF法”。)，进一步优选采用以B3LYP作为泛函数的密度泛函数法(以下称作“B3LYP法”。)。当采用HF法或B3LYP法时，需要指定展开1电子轨道的基函数。当考虑计算时间和计算精度时，可优选采用6-31G (d)基函数、6-31G (d、p)基函数、6-31+G (d、p)基函数。
非结合势能参数进一步分为与范德华势能和库仑势能有关的参数。前者优选以重现密度和蒸发热的实验值的方式来确定。后者的参数是原子所具有的部分电荷，优选以重现由量子化学计算求得的静电势(ESP)的方式来确定。作为重现ESP的方法，可优选采用公知的MK法或CHelpG法。
以上的量子化学计算可以利用可有偿或无偿获得的各种量子化学计算程序包。例如可以优选使用以“Gaussian98(注册商标)”、“ Gaussian03(注册商标)”、“ Gaussian09(注册商标)”、“ GAMESS”这些产品名称市售或公开的通用程序。
当采用分子动力学法与第一原理计算组合的方法作为计算机模拟时，可以对构成原子集合体A的全部原子实施第一原理计算，也可以对构成原子集合体A的一部分原子实施第一原理计算，并且，对于构成原子集合体A的剩余原子可以实施分子力场法。前者不需要与各原子有关的结合信息、势能参数、原子类型，后者只要按照适用第一原理计算的原子范围适当调节与各原子有关的结合信息、势能参数、原子类型等即可。
作为通过计算机模拟产生的系综，例如可以列举：NVE系综(颗粒数N、体积V、势能E一定)、NVT系综(颗粒数N、体积V、温度T一定)、NPT系综(颗粒数N、压力P、体积V一定)等。需要说明的是，系综是指可得到系统的微观状态的集团(统计集团)。
作为边界条件，例如可以列举周期边界条件等。在分子动力学模拟中采用周期边界条件时，在基本单元(Cell)的周围配置图像单元。
温度条件设定为反应式(1)所表示的变化发生时的温度。需要说明的是，温度是指绝对温度，单位是开尔文(Kelvin，K)。
第1坐标获取单元12B获取通过计算机模拟产生的系综中所包含的状态F1～Fi的各状态下的原子集合体A的坐标。状态F1～Fi的各状态是通过实行计算机模拟而产生的微观状态。状态F1～Fi可以是通过计算机模拟而产生的微观状态中的一部分也可以是全部。
i值只要是2以上即可，没有特别限定，可以根据作为计算对象的自由能的种类(例如溶剂化自由能、结合自由能等)适当选择。当计算对象的自由能为溶剂化自由能时，i值优选10000以上，更优选100000以上，进一步优选为1000000以上。当计算对象的自由能为结合自由能时，i值优选100000以上，进一步优选为1000000以上。需要说明的是，对i的上限值没有特别限定，通常为10000000，优选为5000000。
当采用分子动力学模拟作为计算机模拟时，原子集合体A的状态由初期状态(时间T0)随时间而变为状态F1(时间T1)、状态F2(时间T2)、···、状态Fi(时间Ti)。而且，沿着时间序列获取状态F1～Fi的各状态下的原子集合体A的坐标。将所获取的状态F1～Fi (时间T1～Ti)的各状态(各时间)下的原子集合体A的坐标与获取该坐标的各状态(各时间)对应记入，通过控制部12存储在存储部13中。
当采用蒙特卡罗法模拟作为计算机模拟时，通过产生随机数而作出状态F1～Fi，获取状态F1～Fi的各状态下的原子集合体A的坐标。将所获取的状态F1～Fi的各状态下的原子集合体A的坐标与获取该坐标的各状态对应记入，通过控制部12存储在存储部13中。
当采用分子动力学模拟作为计算机模拟时，第1坐标获取单元12B优选在对第1原子集合体模型实行分子动力学模拟之前对其实行能量最小化计算。能量最小化计算例如可以根据存储在存储部13中的第1模拟程序13C的功能来实行。通过能量最小化计算，去除原子集合体模型的初期结构所具有的不自然的结构变形，可以避免计算机模拟的初期阶段的时间积分的发散。
当采用分子动力学模拟作为计算机模拟时，第1坐标获取单元12B优选将第1原子集合体模型(优选为能量最小化计算后的原子集合体模型)平衡化，获取平衡化后的状态F1～Fi下的原子集合体A的坐标。例如，第1坐标获取单元12B对第1原子集合体模型(优选为能量最小化计算后的原子集合体模型)实行计算机模拟，例如在某物理量的值的变动幅度达到恒定的阶段、或者经过了某一定时间的阶段判断为原子集合体模型达到了平衡化，获取平衡化后的状态F1～Fi下的原子集合体A的坐标。
第2坐标获取单元12C根据通过第1坐标获取单元12B获取的原子集合体A的坐标，获取通过片段B与状态F1～Fi的各状态下的原子集合体A结合而产生的原子集合体AB的坐标。
通过第2坐标获取单元12C获取的原子集合体AB的坐标包含构成原子集合体AB的所有原子的坐标。以下，有时将通过片段B与状态Fx(x=1，2，···，i)下的原子集合体A结合而产生的原子集合体AB的坐标记作“RAB(Fx)”。将所获取的原子集合体AB的坐标与作为其基础的原子集合体A的状态对应记入，通过控制部12存储在存储部13中。即，原子集合体AB的坐标RAB(F1)、RAB(F2)、···、RAB(Fi)分别与状态F1、F2、···、Fi对应记入，通过控制部12存储在存储部13中。
通过片段B与状态Fx(x=1、2、···、i)下的原子集合体A结合而产生的原子集合体AB可以是1种原子集合体，也可以是2种以上的原子集合体。即，坐标RAB(Fx)有时是指通过相对于结构a的相对位置不同的第1、第2、···、第p的片段B与状态Fx下的原子集合体A结合而产生的原子集合体AB的坐标。p是从x值中独立选择的2以上的整数。因此，p在x=1、2、···、i中可以是相同的值，在x=1、2、···、i中也可以是不同的值。以下，有时会将第1、第2、···、第p的片段B分别记作“片段B1”、“片段B2”、···、“片段Bp”，而将通过片段B1、B2、···、Bp与状态Fx下的原子集合体A结合而产生的原子集合体AB的坐标分别记作“RAB(Fx，B1)”、“RAB(Fx，B2)”、···、“RAB(Fx，Bp)”。
通过片段B与状态Fx下的原子集合体A结合而产生的原子集合体AB的坐标实际上可以产生原子集合体AB而获取，实际上也可以不产生原子集合体AB而获取。在后一种情况下，例如，确定应该与状态Fx下的原子集合体A结合的片段B相对于结构a的相对位置，根据状态Fx下的原子集合体A的坐标、和所确定的片段B相对于结构a的相对位置，可以获取通过片段B与状态Fx下的原子集合体A结合而产生的原子集合体AB的坐标RAB(Fx)。
对应该与状态Fx下的原子集合体A结合的片段B相对于结构a的相对位置没有特别限定，但从减少具有构成原子集合体A的原子与构成片段B的原子碰撞的结构的原子集合体AB所产生的频率的角度考虑，优选采用相对于构成结构a的原子的坐标限定构成片段B的原子的坐标的约束条件。
约束条件只要能够相对于构成结构a的原子的坐标限定构成片段B的原子的坐标即可，没有特别限定。例如，可以是相对于构成结构a的原子的坐标限定构成片段B的一部分原子的坐标的条件，也可以是相对于构成结构a的原子的坐标限定构成片段B的所有原子的坐标的条件。
在一实施方式中，第2坐标获取单元12C以x=1～i反复实行以下的处理：
处理W101，关于应该与状态Fx下的原子集合体A结合的片段B，规定选自构成片段B的原子的原子b相对于结构a的相对位置条件(例如相对于构成结构a的1个或2个以上的原子的距离、角度等)；
处理W102，根据所规定的相对位置条件，获取原子b的坐标；
处理W103，根据所获取的原子b的坐标和状态Fx下的原子集合体A的坐标RA(Fx)，获取通过原子b与状态Fx下的原子集合体A结合而产生的原子集合体Ab的坐标RAb(Fx)；以及
处理W104，以原子b的坐标为基准，通过计算机模拟、公知的构象产生方法等，产生构成片段B的原子b以外的原子的坐标，获取通过片段B与状态Fx下的原子集合体A结合而产生的原子集合体AB的坐标RAB(Fx)。由此，可以获取通过限定了相对于结构a的相对位置的片段B与状态F1～Fi的各状态下的原子集合体A结合而产生的原子集合体AB的坐标RAB(F1)～RAB(Fi)。
在处理W101中，关于应该与状态Fx下的原子集合体A结合的片段B，第2坐标获取单元12C可以规定彼此不同的第1、第2、···、第p的相对位置条件。p是从x的值中独立选择的2以上的整数。因此，在x=1、2、···、i中p可以是相同的值，在x=1、2、···、i中p也可以是不同的值。
在处理W101中，关于应该与状态Fx下的原子集合体A结合的片段B，当规定彼此不同的第1、第2、···、第p的相对位置条件时，关于第1～第p的各相对位置条件，第2坐标获取单元12C反复进行处理W102～W104。由此，可以获取通过相对于结构a的相对位置不同的片段B1、B2、···、Bp与状态Fx下的原子集合体A结合而产生的原子集合体AB的坐标RAB(Fx，B1)、RAB(Fx，B2)、···、RAB(Fx，Bp)。
在另一实施方式中，第2坐标获取单元12C实行以下的处理：
处理W201，制作将由结构a和与该结构a结合的片段B构成或者包含结构a和与该结构a结合的片段B而构成的原子集合体C进行模型化而获得的第3原子集合体模型；以及
处理W202，通过针对所制作的第3原子集合体模型的计算机模拟，获取第1～第k的状态H1～Hk (k为2以上的整数。)的各状态下的原子集合体C的坐标(以下，有时将状态H1、H2、···、Hk下的原子集合体C的坐标分别记作“RC(H1)”、“ RC(H2)”、···、“ RC(Hk)”。)，
之后，以x=1～i反复实行以下的处理：
处理W203，选择由选自构成结构a的原子的1个或2个以上的原子构成的选择原子组(需要说明的是，不仅包含选择2个以上的原子的情形，还包含选择1个原子的情形，为方便起见称作“选择原子组”。)，使选自状态H1～Hk的1种或2种以上的状态下的原子集合体C的选择原子组的坐标相对于状态Fx下的原子集合体A的选择原子组的坐标旋转和/或并进，从而在原子集合体A的选择原子组与原子集合体C的选择原子组之间制作对应的原子间距离的平方和达到最小的原子集合体C的坐标(使原子集合体C的选择原子组的坐标与原子集合体A的选择原子组的坐标拟合)，根据所制作的原子集合体C的坐标，使选自状态H1～Hk的1种或2种以上的状态下的原子集合体C与原子集合体A重叠；以及
处理W204，根据状态Fx下的原子集合体A的坐标RA(Fx)、以及与状态Fx下的原子集合体A重叠的原子集合体C的片段B的1个或2个以上的坐标，获取通过片段B与状态Fx下的原子集合体A结合而产生的原子集合体AB的坐标。由此，可以获取通过限定了相对于结构a的相对位置的片段B与状态F1～Fi的各状态下的原子集合体A结合而产生的原子集合体AB的坐标RAB(F1)～RAB(Fi)。
原子集合体C可以仅由分子C’构成，所述分子C’由结构a和与该结构a结合的片段B构成(此时分子C’的周围为真空)，也可以由分子C’和1种或2种以上的其它分子(例如水分子等溶剂分子)构成，所述分子C’由结构a和与该结构a结合的片段B构成。
在处理W201中，第2坐标获取单元12C制作将原子集合体C进行模型化而获得的第3原子集合体模型。通过第2坐标获取单元12C制作的第3原子集合体模型的相关数据(例如构成原子集合体C的各原子的坐标、种类、质量、部分电荷、原子间的结合信息等)通过控制部12存储在存储部13中。
第2坐标获取单元12C在制作第3原子集合体模型时使用存储在存储部13中的原子集合体C的相关数据(未图示)。原子集合体C的相关数据例如以控制部12可读取的文件形式存储在存储部13中。原子集合体C的相关数据只要能够制作将原子集合体C进行模型化而获得第3原子集合体模型即可，没有特别限定。作为原子集合体C的相关数据，例如可以列举：分子C’的相关数据(例如构成分子C’的各原子的坐标、种类、质量、部分电荷、原子间的结合信息等)、分子C’的周边环境的相关数据(例如位于分子C’周围的溶剂分子的有无、种类、个数、构成溶剂分子的原子的坐标等)等。
第3原子集合体模型只要能够进行计算机模拟即可，没有特别限定。作为原子集合体模型，可以列举：全原子模型、珠-簧模型、联合原子模型等。
针对第3原子集合体模型的计算机模拟只要基于统计力学理论即可，没有特别限定。作为代表性的计算机模拟，例如可以列举：分子动力学法、蒙特卡罗法、将分子动力学法或蒙特卡罗法与第一原理计算组合的方法等。针对第3原子集合体模型的计算机模拟可以按照与针对第1原子集合体模型的计算机模拟相同的方式实行。
当采用分子动力学模拟作为计算机模拟时，第2坐标获取单元12C优选在对第3原子集合体模型实行分子动力学模拟之前对其实行能量最小化计算。能量最小化计算例如可以根据存储在存储部13中的第1模拟程序13C的功能来实行。通过能量最小化计算，去除原子集合体模型的初期结构所具有的不自然的结构变形，可以避免计算机模拟的初期阶段的时间积分的发散。
当采用分子动力学法作为计算机模拟时，第2坐标获取单元12C优选将第3原子集合体模型(优选为能量最小化计算后的原子集合体模型)平衡化，获取平衡化后的状态H1～Hk下的原子集合体C的坐标。例如，第2坐标获取单元12C对第3原子集合体模型(优选为能量最小化计算后的原子集合体模型)实行计算机模拟，例如在某物理量的值的变动幅度达到恒定的阶段、或者经过了某一定时间的阶段判断为原子集合体模型达到了平衡化，获取平衡化后的状态H1～Hk下的原子集合体C的坐标。
在处理W203中，相对于状态Fx下的原子集合体A的坐标RA(Fx)，第2坐标获取单元12C可以从状态H1～Hk下的原子集合体C的坐标中选择片段B相对于结构a的相对位置彼此不同的p组的坐标。p是从x的值中独立选择的2以上的整数。因此，在x=1、2、···、i中p可以是相同的值，在x=1、2、···、i中p也可以是不同的值。
在处理W203中，从状态H1～Hk下的原子集合体C的坐标中选择片段B相对于结构a的相对位置彼此不同的p组的坐标时，第2坐标获取单元12C关于p组的各坐标反复进行处理W203～W204。由此，可以获取通过相对于结构a的相对位置不同的片段B1、B2、···、Bp与状态Fx下的原子集合体A结合而产生的原子集合体AB的坐标RAB(Fx，B1)、RAB(Fx，B2)、···、RAB(Fx，Bp)。
在处理W203中，以选择原子组的坐标为基准，使选自状态H1～Hk的1种或2种以上状态下的原子集合体C的坐标相对于状态Fx下的原子集合体A的坐标RA(Fx)旋转和/或并进，从而使选自状态H1～Hk的1种或2种以上状态下的原子集合体C与原子集合体A重叠，上述方法优选最小平方法。通过采用该最小平方法，在选自状态H1～Hk的1种或2种以上的状态下的原子集合体C的坐标中，可以维持构成选择原子组的原子与构成片段B的原子的相对位置关系，同时在选自状态H1～Hk的1种或2种以上的状态下的原子集合体C的坐标中，可以向状态Fx下的原子集合体A的坐标RA(Fx)中添加构成片段B的原子的坐标。由此，可以减少构成片段B的原子与构成结构a的原子碰撞的结构的产生频率，可以使片段B与原子集合体A结合。即，可以提高抽样效率，在统计学上可以以高精度获得计算结果。
对最小平方法的具体实施方法没有特别限定。例如，可以如下实施最小平方法。在状态F1～Fi的各状态下的原子集合体A的坐标中将选择原子组的坐标记作“(xs_A：n、ys_A：n、zs_A：n)”、而在选自状态H1～Hk的1种状态下的原子集合体C的坐标中将选择原子组的坐标记作“(xC：n、yC：n、zC：n)”时，在原子集合体C的坐标中保持选择原子组的原子间的相对坐标(内部坐标)的条件下，使原子集合体C的选择原子组的坐标旋转和/或并进，从而求出(xT：n、yT：n、zT：n)，使数学式(7)达到最小。这里，以坐标的附带文字形式赋予的n是对构成选择原子组的各个原子赋予的数字，坐标(xs_A：n、ys_A：n、zs_A：n)和(xT：n、yT：n、zT：n)表示构成选择原子组的各个原子的坐标。而且，数学式(7)的Cn为零以上的实数，可以对构成原子集合体A和原子集合体C中的选择原子组的各原子进行设定。根据对构成选择原子组的各原子设定的Cn的相对比，可以设定通过基于数学式(7)的最小平方法获得的构成选择原子组的原子的坐标一致的重要度。例如，当构成选择原子组的原子包含重原子(氢原子以外的原子)和氢原子两者时，与氢原子相比，将该重原子的Cn设定得大，从而可以使原子集合体A和原子集合体C中的选择原子组的该重原子的坐标重点一致。
[数学式9]
接下来，在通过最小平方法获得的坐标(xT：n、yT：n、zT：n)中，以构成原子集合体C中的选择原子组的原子的坐标一致的方式，使构成原子集合体C的原子旋转和/或并进，从而保持构成原子集合体C的原子间的相对坐标(内部坐标)，同时使构成原子集合体C的原子的坐标移动，将构成该原子集合体C中的片段B的原子的坐标(其中，点电荷等虚拟原子除外)作为构成原子集合体AB中的片段B的原子的坐标，从而可以获取原子集合体AB的坐标RAB(F1)～RAB(Fi)。
当片段B包含点电荷(虚拟原子)时，第2坐标获取单元12C优选在构成原子集合体A的结构a的原子的电荷参数中添加片段B的点电荷。在计算机模拟上，使构成除虚拟原子以外的片段B的原子与构成原子集合体A的结构a的原子结合时，有时需要改变构成原子集合体A的结构a的原子的部分电荷(电荷参数)。例如，当构成去除了虚拟原子的片段B的原子构成吸电子基团或供电子基团、并使构成片段B的原子与构成原子集合体A的结构a的原子结合时，在构成原子集合体A的结构a的原子与构成去除了虚拟原子的片段B的原子之间生成键的情况下，通常在构成原子集合体A的结构a的原子中，与片段B结合的原子和/或与片段B结合的原子周围的原子(例如，在构成原子集合体A的原子中，与构成去除了虚拟原子的片段B的原子结合的原子所结合的原子)的部分电荷(电荷参数)会发生变化。还考虑了将构成部分电荷(电荷参数)因是否存在与片段B的结合而发生变化的结构a的原子作为构成片段B的原子(实际存在原子)来处理、而不作为构成结构a的原子来处理的方法，但通过包含具有与该原子的部分电荷(电荷参数)的变化量同等的电荷的点电荷作为构成片段B的原子(虚拟原子)，可以减小原子集合体AB(终状态)与原子集合体A(始状态)的变化，因此可以提高计算精度。
第1相互作用能φ频率分布制作单元12D根据通过第2坐标获取单元12C获取的原子集合体AB的坐标，计算结构a和与该结构a结合的片段B的相互作用能φ，制作表示相互作用能φ的各级频率的频率分布。所计算的相互作用能φ和所制作的频率分布通过控制部12存储在存储部13中。将相互作用能φ的各级与作为其计算基础的原子集合体AB的坐标对应记入，通过控制部12存储在存储部13中。
相互作用能φ与原子集合体AB的各坐标RAB(F1)～RAB(Fi)相关联进行计算。相互作用能φ例如可以根据存储在存储部13中的第1模拟程序13C的功能来计算。在计算相互作用能φ时，对使用的模拟程序和势能参数没有特别限定，但由于原子集合体AB的坐标是通过第1坐标获取单元12B所实行的计算机模拟而求得的，因此优选使用与第1坐标获取单元12B所实行的计算机模拟中使用的模拟程序和势能参数相同者。
表示相互作用能φ的各级频率的频率分布例如可以以直方图的形式来制作，其中横轴表示相互作用能φ的各级、而纵轴表示相互作用能φ的各级频率。
相互作用能φ的各级例如可以以相互作用能φ为中心、以具有相互作用能间隔Δφ的相互作用能区间[φ-Δφ/2～φ+Δφ/2]作为级间隔来制作，对于所有区间而言，各相互作用能区间中的Δφ均可以是一定的Δφ，也可以使用根据相互作用能区间而适当变更了的Δφ。在进行数学式(1)的第2项的计算上，优选使P0(φ)的级间隔Δφ与P(φ)的级间隔Δφ一致，并优选使P(φ)的级间隔Δφ与P0(φ)的级间隔Δφ一致。对相互作用能φ的分配函数(分割数，partition function)、即级间隔Δφ的个数没有特别限定，但优选50～500、进一步优选为250～500。
在数学式(1)的第2项中，包含用相互作用能φ的出现概率P0(φ)除以相互作用能φ的出现概率P(φ)的计算过程，因此P(φ)是不同于零的值，若存在P0(φ)为零的相互作用能区间，则数学式(1)的第2项发散，无法根据数学式(1)来计算自由能。因此，在所有的相互作用能区间[φ-Δφ/2～φ+Δφ/2]，优选在每个相互作用能区间适当设定Δφ，使P0(φ)达到一定。通过在每个相互作用能区间适当设定Δφ、并使P(φ)的级间隔Δφ与P0(φ)的级间隔Δφ一致，可以在统计学上获得精度高的计算结果。
第1相互作用能φ出现概率计算单元12E根据通过第1相互作用能φ频率分布制作单元12D制作的频率分布，计算相互作用能φ的各级的出现概率P0(φ)。所计算的P0(φ)通过控制部12存储在存储部13中。
第1相互作用能φ出现概率计算单元12E通过将表示相互作用能φ的各级频率的频率分布标准化，可以计算相互作用能φ的各级的出现概率P0(φ)。这里，标准化是指用相互作用能φ的各级频率总和除以频率分布中的相互作用能φ的各级频率。
第1相互作用能ε频率分布制作单元12F根据通过第2坐标获取单元12C获取的原子集合体AB的坐标，计算通过从原子集合体AB中去除由结构a和与该结构a结合的片段B构成的结构aB而产生的原子集合体的一部分或全部与片段B的相互作用能ε，制作通过第1相互作用能φ频率分布制作单元12D制作的频率分布的相互作用能φ的各级中的、表示相互作用能ε的各级频率的频率分布。所计算的相互作用能ε和所制作的频率分布通过控制部12存储在存储部13中。
在一实施方式中，第1相互作用能ε频率分布制作单元12F根据通过第2坐标获取单元12C获取的原子集合体AB的坐标，关于原子集合体AB的各坐标RAB(F1)～RAB(Fi)，计算通过从原子集合体AB中去除由结构a和与该结构a结合的片段B构成的结构aB而产生的原子集合体的一部分或全部与片段B的相互作用能ε，之后，根据通过第1相互作用能φ频率分布制作单元12D制作的频率分布，抽取相互作用能φ的各级中的、通过从原子集合体AB中去除由结构a和与该结构a结合的片段B构成的结构aB而产生的原子集合体的一部分或全部与片段B的相互作用能ε，制作通过第1相互作用能φ频率分布制作单元12D制作的频率分布的相互作用能φ的各级中的表示相互作用能ε的各级频率的频率分布。
在又一实施方式中，第1相互作用能ε频率分布制作单元12F根据通过第1相互作用能φ频率分布制作单元12D制作的频率分布，从通过第2坐标获取单元12C获取的原子集合体AB的坐标中抽取属于相互作用能φ的各级的原子集合体AB的坐标，之后，根据所抽取的原子集合体AB的坐标，计算相互作用能φ的各级中的、通过从原子集合体AB中去除由结构a和与该结构a结合的片段B构成的结构aB而产生的原子集合体的一部分或全部与片段B的相互作用能ε，制作通过第1相互作用能φ频率分布制作单元12D制作的频率分布的相互作用能φ的各级中的表示相互作用能ε的各级频率的频率分布。
在原子集合体AB由配体和溶剂分子(例如水分子)构成、且由结构a和与该结构a结合的片段B构成的结构aB为该配体的情况下，通过从原子集合体AB中去除由结构a和与该结构a结合的片段B构成的结构aB而产生的原子集合体的一部分或全部与片段B的相互作用能ε是指片段B与存在于该片段B周边的各水分子的相互作用能。另外，在原子集合体AB由配体、蛋白和水分子构成、且由结构a和与该结构a结合的片段B构成的结构aB为该配体的情况下，通过从原子集合体AB中去除由结构a和与该结构a结合的片段B构成的结构aB而产生的原子集合体的一部分或全部与片段B的相互作用能ε是指片段B与各水分子间的相互作用能和片段B与蛋白间的相互作用能。这里列举的相互作用能ε的具体例子不仅适用于通过第1相互作用能ε频率分布制作单元12F计算的相互作用能ε，对于通过第2相互作用能ε频率分布制作单元12L计算的相互作用能ε也适用。
相互作用能ε与通过第1相互作用能φ频率分布制作单元12D制作的频率分布(例如直方图)中所含的相互作用能φ的所有级相关联进行计算。相互作用能ε例如可以根据存储在存储部13中的第1模拟程序13C的功能来计算。在计算相互作用能ε时，对使用的模拟程序和势能参数没有特别限定，但由于原子集合体AB的坐标是通过第1坐标获取单元12B所实行的计算机模拟而求得的，因此优选使用与第1坐标获取单元12B所实行的计算机模拟中使用的模拟程序和势能参数相同者。
表示相互作用能ε的各级频率的频率分布与相互作用能φ的各级中的、相互作用能ε相关联进行制作。表示相互作用能ε的各级频率的频率分布例如可以以直方图的形式制作，其中横轴表示相互作用能ε的各级、纵轴表示相互作用能ε的各级频率。
相互作用能ε的各级例如可以以相互作用能ε为中心、以具有相互作用能间隔Δε的相互作用能区间[ε-Δε/2～ε+Δε/2]作为级间隔而制作，对于所有区间而言，各相互作用能区间中的Δε均可以是一定的Δε，也可以使用根据相互作用能区间而适当变更了的Δε。因相互作用能ε的分配函数、即级间隔Δε的个数取决于通过模拟得到的抽点打印数目、片段B周围的水分子等的数目、相互作用能ε的能量区间宽度等，故没有特别限定，但优选500～5000，进一步优选为2000～5000。
第1相互作用能ε出现概率计算单元12G根据通过第1相互作用能ε频率分布制作单元12F制作的频率分布，计算相互作用能φ的各级中的相互作用能ε的各级的出现概率P0’(ε；φ)。所计算的P0’(ε；φ)通过控制部12存储在存储部13中。需要说明的是，P0’(ε；φ)中的“φ”是用于明确在相互作用能φ的每个等级均计算相互作用能ε的各级的出现概率P0’(ε；φ)的符号。
第1相互作用能ε出现概率计算单元12G通过将表示相互作用能ε的各级频率的频率分布标准化，可以计算相互作用能ε的各级的出现概率P0’(ε；φ)。这里，标准化是指用相互作用能ε的各级频率的总和除以频率分布中的相互作用能ε的各级频率。
第2原子集合体模型制作单元12H制作将变化后的原子集合体AB进行模型化而获得的第2原子集合体模型。通过第2原子集合体模型制作单元12H制作的第2原子集合体模型的相关数据(例如构成原子集合体AB的各原子的坐标、种类、质量、部分电荷、原子间的结合信息等)通过控制部12存储在存储部13中。
第2原子集合体模型制作单元12H在制作第2原子集合体模型时，使用存储在存储部13中的变化后的原子集合体AB的相关数据13D。变化后的原子集合体AB的相关数据13D例如以控制部12可读取的文件形式存储在存储部13中。变化后的原子集合体AB的相关数据13D只要能够制作将变化后的原子集合体AB进行模型化而获得的第2原子集合体模型即可，没有特别限定。作为变化后的原子集合体AB的相关数据13D，例如可以列举：构成原子集合体AB的各原子的坐标、种类、质量、部分电荷、原子间的结合信息等。
第2原子集合体模型只要可以进行计算机模拟即可，没有特别限定。作为原子集合体模型，可以列举：全原子模型、珠-簧模型、联合原子模型等。需要说明的是，珠-簧模型是指将构成原子集合体AB的单体单元视为1个珠粒(段)、并通过假想的弹簧使之结合的模型，而联合原子模型是指将氢原子包含在重原子(例如碳原子)中作为1个原子(质点)进行处理的模型。
第3坐标获取单元12I通过以针对通过第2原子集合体模型制作单元12H制作的第2原子集合体模型的计算机模拟的结果的形式输出的抽点打印，获取第1～第j的状态G1～Gj (j为2以上的整数。)的各状态下的原子集合体AB的坐标。这里，抽点打印包含构成原子集合体AB的所有原子的坐标。即，在状态G1～Gj的各状态下的抽点打印中，分别包含构成状态G1～Gj的各状态下的原子集合体AB的所有原子的坐标。以下，有时将状态G1、G2、···、Gi下的原子集合体AB的坐标分别称作坐标RAB(G1)、RAB(G2)、···、RAB(Gj)。通过第3坐标获取单元获取的原子集合体AB的坐标通过控制部12存储在存储部13中。
第3坐标获取单元12I所实行的计算机模拟只要基于统计力学理论即可，没有特别限定。作为代表性的计算机模拟，例如可以列举：分子动力学法、蒙特卡罗法、将分子动力学法或蒙特卡罗法与第一原理计算组合的方法等。
第3坐标获取单元12I根据存储在存储部13中的第2原子集合体模型的相关数据、第2模拟条件的相关数据13E、第2模拟程序13F等实行计算机模拟。模拟条件的相关数据13E以控制部12可读取的文件形式存储在存储部13中。
模拟条件只要能够实行计算机模拟即可，没有特别限定。关于计算机模拟的具体说明与第1坐标获取单元12B所实行的计算机模拟相同，因此省略。
第3坐标获取单元12I获取通过计算机模拟而产生的系综所包含的状态G1～Gj的各状态下的原子集合体AB的坐标。状态G1～Gj的各状态是指通过实行计算机模拟而产生的微观状态。状态G1～Gj可以是通过计算机模拟而产生的微观状态中的一部分，也可以是其全部。
j的值只要是2以上即可，没有特别限定，可以根据作为计算对象的自由能的种类(例如溶剂化自由能、结合自由能等)而适当选择。当计算对象的自由能为溶剂化自由能时，j的值优选10000以上，更优选为100000以上。当计算对象的自由能为结合自由能时，j的值优选100000以上，进一步优选为1000000以上。需要说明的是，对j的上限值没有特别限定，但通常为10000000，优选为5000000。
当采用分子动力学法作为计算机模拟时，原子集合体AB的状态由初期状态(时间T0)随时间而变为状态G1(时间T1)、状态G2(时间T2)、···、状态Gj(时间Tj)。而且，沿时间序列获取状态G1～Gj的各状态下的原子集合体AB的坐标。将所获取的状态G1～Gj (时间T1～Tj)的各状态(各时间)下的原子集合体AB的坐标与获取该坐标的各状态(各时间)对应记入，通过控制部12存储在存储部13中。
当采用蒙特卡罗法作为计算机模拟时，通过产生随机数而作出状态G1～Gj，获取状态G1～Gj的各状态下的原子集合体AB的坐标。将所获取的状态G1～Gj的各状态下的原子集合体AB的坐标与获取该坐标的各状态对应记入，通过控制部12存储在存储部13中。
当采用分子动力学法作为计算机模拟时，第3坐标获取单元12I优选在对第2原子集合体模型实行分子动力学模拟之前对其实行能量最小化计算。能量最小化计算例如可以根据存储在存储部13中的第2模拟程序13F的功能来实行。通过能量最小化计算，去除原子集合体模型的初期结构所具有的不自然的结构变形，可以避免计算机模拟的初期阶段的时间积分的发散。
当采用分子动力学法作为计算机模拟时，第3坐标获取单元12I优选将第2原子集合体模型(优选能量最小化计算后的原子集合体模型)进行平衡化，获取平衡化后的状态G1～Gj下的原子集合体AB的坐标。例如，第3坐标获取单元12I对第2原子集合体模型(优选为能量最小化计算后的原子集合体模型)实行计算机模拟，例如在某物理量达到了阈值的阶段、或者在经过了某一定时间的阶段判断为原子集合体模型达到了平衡化，获取平衡化后的状态G1～Gj下的原子集合体AB的坐标。
第2相互作用能φ频率分布制作单元12J根据通过第3坐标获取单元12I获取的原子集合体AB的坐标，计算结构a和与该结构a结合的片段B的相互作用能φ，制作表示相互作用能φ的各级频率的频率分布。将相互作用能φ的各级与作为其计算基础的原子集合体AB的坐标对应记入，通过控制部12存储在存储部13中。
相互作用能φ与原子集合体AB的各坐标RAB(G1)～RAB(Gj)相关联进行计算。相互作用能φ例如可以根据第2模拟程序13F的功能来计算。在计算相互作用能φ时，对使用的模拟程序和势能参数没有特别限定，但由于原子集合体AB的坐标是通过第3坐标获取单元12I所实行的计算机模拟而求得的，因此优选使用与第3坐标获取单元12I所实行的计算机模拟中使用的模拟程序和势能参数相同者。
表示相互作用能φ的各级频率的频率分布例如可以以直方图的形式制作，其中横轴表示相互作用能φ的各级、纵轴表示相互作用能φ的各级频率。
相互作用能φ的各级例如可以以相互作用能φ为中心、以具有相互作用能间隔Δφ的相互作用能区间[φ-Δφ/2～φ+Δφ/2]作为级间隔来制作，各相互作用能区间中的Δφ可以是对于所有的区间而言均为一定的Δφ，也可以使用根据相互作用能区间而适当变更了的Δφ。在进行数学式(1)的第2项的计算上，优选使P0(φ)的级间隔Δφ与P(φ)的级间隔Δφ一致，进一步优选使P(φ)的级间隔Δφ与P0(φ)的级间隔Δφ一致。对相互作用能φ的分配函数、即级间隔Δφ的个数没有特别限定，但优选50～500、进一步优选为250～500。
在数学式(1)的第2项中，由于包含用相互作用能出现概率P0(φ)除以相互作用能存在概率P(φ)的计算过程，所以P(φ)为不同于零的值，若存在相互作用能出现概率P0(φ)为零的相互作用能区间，则数学式(1)的第2项发散，无法根据数学式(1)来计算自由能。因此，优选在所有的相互作用能区间[φ-Δφ/2～φ+Δφ/2]中，在每个相互作用能区间适当设定Δφ，使P0(φ)达到一定。通过在每个相互作用能区间适当设定Δφ、并使P(φ)的级间隔Δφ与P0(φ)的级间隔Δφ一致，可以在统计学上得到精度高的计算结果。
第2相互作用能φ出现概率计算单元12K根据通过第2相互作用能φ频率分布制作单元12J制作的频率分布，计算相互作用能φ的各级的出现概率P(φ)。所计算的P(φ)通过控制部12存储在存储部13中。
第2相互作用能φ出现概率计算单元12K通过将表示相互作用能φ的各级频率的频率分布标准化，可以计算相互作用能φ的各级的出现概率P(φ)。这里，标准化是指用相互作用能φ的各级频率的总和除以频率分布中的相互作用能φ的各级频率。
第2相互作用能ε频率分布制作单元12L根据通过第3坐标获取单元12I获取的原子集合体AB的坐标，计算通过从原子集合体AB中去除由结构a和与该结构a结合的片段B构成的结构aB而产生的原子集合体的一部分或全部与片段B的相互作用能ε，制作表示相互作用能ε的各级频率的频率分布。所计算的相互作用能ε和所制作的频率分布通过控制部12存储在存储部13中。
第2相互作用能ε频率分布制作单元12L可以将相互作用能ε与通过第2相互作用能φ频率分布制作单元12J制作的频率分布中的相互作用能φ的各级无关联地进行计算，也可以关联着进行计算。将相互作用能ε与相互作用能φ的各级关联着进行计算时，相互作用能ε与通过第2相互作用能φ频率分布制作单元12J制作的频率分布(例如直方图)中所含的相互作用能φ的所有级相关联进行计算，表示相互作用能ε的各级频率的频率分布与相互作用能φ的各级中的相互作用能ε相关联进行制作。
在将相互作用能ε与相互作用能φ的各级无关联地进行计算时，第2相互作用能ε频率分布制作单元12L根据通过第3坐标获取单元12I获取的原子集合体AB的坐标，关于原子集合体AB的各坐标RAB(G1)～RAB(Gj)，计算通过从原子集合体AB中去除由结构a和与该结构a结合的片段B构成的结构aB而产生的原子集合体的一部分或全部与片段B的相互作用能ε，制作表示相互作用能ε的各级频率的频率分布。
在将相互作用能ε与相互作用能φ的各级关联着进行计算时，第2相互作用能ε频率分布制作单元12L根据通过第3坐标获取单元12I获取的原子集合体AB的坐标，计算通过从原子集合体AB中去除由结构a和与该结构a结合的片段B构成的结构aB而产生的原子集合体的一部分或全部与片段B的相互作用能ε，制作通过第2相互作用能φ频率分布制作单元12J制作的频率分布的相互作用能φ的各级中的、表示相互作用能ε的各级频率的频率分布。
在一实施方式中，第2相互作用能ε频率分布制作单元12L根据通过第3坐标获取单元12I获取的原子集合体AB的坐标，关于原子集合体AB的各坐标RAB(G1)～RAB(Gj)，计算通过从原子集合体AB中去除由结构a和与该结构a结合的片段B构成的结构aB而产生的原子集合体的一部分或全部与片段B的相互作用能ε，之后，根据通过第2相互作用能φ频率分布制作单元12J制作的频率分布，抽取相互作用能φ的各级中的、通过从原子集合体AB中去除由结构a和与该结构a结合的片段B构成的结构aB而产生的原子集合体的一部分或全部与片段B的相互作用能ε，制作通过第2相互作用能φ频率分布制作单元12J制作的频率分布的相互作用能φ的各级中的、表示相互作用能ε的各级频率的频率分布。
在又一实施方式中，第2相互作用能ε频率分布制作单元12L根据通过第2相互作用能φ频率分布制作单元12J制作的频率分布，从通过第3坐标获取单元12I获取的原子集合体AB的坐标中抽取属于相互作用能φ的各级的原子集合体AB的坐标，之后，根据所抽取的原子集合体AB的坐标，计算相互作用能φ的各级中的、通过从原子集合体AB中去除由结构a和与该结构a结合的片段B构成的结构aB而产生的原子集合体的一部分或全部与片段B的相互作用能ε，制作通过第2相互作用能φ频率分布制作单元12J制作的频率分布的相互作用能φ的各级中的、表示相互作用能ε的各级频率的频率分布。
相互作用能ε例如可以根据第2模拟程序13F的功能来计算。在计算相互作用能ε时，对使用的模拟程序和势能参数没有特别限定，但由于原子集合体AB的坐标是通过第3坐标获取单元12I所实行的计算机模拟而求得的，所以优选使用与第3坐标获取单元12I所实行的计算机模拟中使用的模拟程序和势能参数相同者。
表示相互作用能ε的各级频率的频率分布例如可以以直方图的形式制作，其中横轴表示相互作用能ε的各级、纵轴表示相互作用能ε的各级频率。
相互作用能ε的各级例如可以以相互作用能ε为中心、以具有相互作用能间隔Δε的相互作用能区间[ε-Δε/2～ε+Δε/2]作为级间隔来制作，各相互作用能区间中的Δε可以是对于所有的区间而言均为一定的Δε，也可以使用根据相互作用能区间而适当变更了的Δε。另外，P0’(ε；φ)中的级间隔Δε与P’(ε)中的级间隔Δε优选一致。由于相互作用能ε的分配函数、即级间隔Δε的个数取决于通过模拟而得到的抽点打印数目、片段B周围的水分子的数目、相互作用能ε的能量区间宽度等，因此没有特别限定，但优选500～5000、进一步优选为2000～5000。
第2相互作用能ε出现概率计算单元12M根据通过第2相互作用能ε频率分布制作单元12L制作的频率分布，计算相互作用能ε的各级的出现概率P’(ε)。所计算的P’(ε)通过控制部12存储在存储部13中。
第2相互作用能ε出现概率计算单元12M通过将表示相互作用能ε的各级频率的频率分布标准化，可以计算相互作用能ε的各级的出现概率P’(ε)。这里，标准化是指用相互作用能ε的各级频率的总和除以频率分布中的相互作用能ε的各级频率。
第2相互作用能ε频率分布制作单元12L将相互作用能ε与通过第2相互作用能φ频率分布制作单元12J制作的频率分布中的相互作用能φ的各级无关联地进行计算时，第2相互作用能ε出现概率计算单元12M根据通过第2相互作用能ε频率分布制作单元12L制作的频率分布，而与相互作用能φ的各级无关联，计算相互作用能ε的各级的出现概率P’(ε)。
当第2相互作用能ε频率分布制作单元12L将相互作用能ε与通过第2相互作用能φ频率分布制作单元12J制作的频率分布中的相互作用能φ的各级关联着进行计算时，第2相互作用能ε出现概率计算单元12M根据通过第2相互作用能ε频率分布制作单元12L制作的频率分布，计算相互作用能φ的各级中的相互作用能ε的各级的出现概率P’(ε；φ)。需要说明的是，P’(ε；φ)中的“φ”是用于明确在相互作用能φ的每个等级计算相互作用能ε的各级的出现概率P’(ε；φ)的符号。
∫Δν(φ)P(φ)dφ计算单元12N根据通过第2相互作用能φ出现概率计算单元12K计算的P(φ)、通过第1相互作用能ε出现概率计算单元12G计算的P0’(ε；φ)、以及通过第2相互作用能ε出现概率计算单元12M计算的P’(ε)，计算相互作用能ε所引起的自由能变化量∫Δν(φ)P(φ)dφ。
“∫Δν(φ)P(φ)dφ”中的“Δν(φ)”表示相互作用能φ的各级中的相互作用能ε所引起的自由能变化量。当原子集合体AB由配体和溶剂分子(例如水分子)构成、且由结构a和与该结构a结合的片段B构成的结构aB为该配体时，相互作用能ε所引起的自由能变化量是指由片段B与存在于该片段B周边的水分子的相互作用能所引起的自由能变化量。另外，当原子集合体AB由配体、蛋白和水分子构成、且由结构a和与该结构a结合的片段B构成的结构aB为该配体时，相互作用能ε所引起的自由能变化量是指片段B与水分子间的相互作用能所引起的自由能变化量和片段B与蛋白间的相互作用能所引起的自由能变化量的和。
∫Δν(φ)P(φ)dφ计算单元12N优选利用能量表示法，根据通过第2相互作用能φ出现概率计算单元12K计算的P(φ)、通过第1相互作用能ε出现概率计算单元12G计算的P0’(ε；φ)、以及通过第2相互作用能ε出现概率计算单元12M计算的P’(ε)，计算相互作用能ε所引起的自由能变化量∫Δν(φ)P(φ)dφ。
在一实施方式中，∫Δν(φ)P(φ)dφ计算单元12N不计算Δν(φ)而根据通过第2相互作用能φ出现概率计算单元12K计算的P(φ)、通过第1相互作用能ε出现概率计算单元12G计算的P0’(ε；φ)、以及通过第2相互作用能ε出现概率计算单元12M计算的P’(ε)，利用能量表示法计算∫Δν(φ)P(φ)dφ。这里，使用的P’(ε)是指与相互作用能φ的各级无关联地进行计算的相互作用能ε的各级的出现概率。
在又一实施方式中，∫Δν(φ)P(φ)dφ计算单元12N根据通过第1相互作用能ε出现概率计算单元12G计算的P0’(ε；φ)、以及通过第2相互作用能ε出现概率计算单元12M计算的P’(ε；φ)，计算相互作用能φ的各级中的、相互作用能ε所引起的自由能变化量Δν(φ)，根据所计算的Δν(φ)和通过第2相互作用能φ出现概率计算单元12K计算的P(φ)来计算∫Δν(φ)P(φ)dφ。此时，∫Δν(φ)P(φ)dφ计算单元12N优选利用能量表示法，根据通过第1相互作用能ε出现概率计算单元12G计算的P0’(ε；φ)、以及通过第2相互作用能ε出现概率计算单元12M计算的P’(ε；φ)，计算相互作用能φ的各级中的、相互作用能ε所引起的自由能变化量Δν(φ)。
ΔG计算单元12O根据通过第1相互作用能φ出现概率计算单元12E计算的P0(φ)、通过第2相互作用能φ出现概率计算单元12K计算的P(φ)、通过∫Δν(φ)P(φ)dφ计算单元12N计算的∫Δν(φ)P(φ)dφ和数学式(1)来计算ΔG，
[数学式10]
[式中，R表示气体常数，T表示反应式(1)所表示的变化发生的绝对温度。]。
以下，根据附图，对本发明的一实施方式所涉及的计算处理次序进行说明。
图2是显示本实施方式所涉及的计算装置1的处理次序的流程图。
控制部12在实行变化前的原子集合体A的相关计算处理(步骤S110～190)和变化后的原子集合体AB的相关计算处理(步骤S210～280)时，可以在实行其中一种处理后再实行另一种处理，也可以并行实行这两种处理。
以下，对步骤S110～S150进行说明。
控制部12依次实施步骤S110～S150。
在步骤S110中，控制部12利用存储在存储部13中的变化前的原子集合体A的相关数据13A，制作将变化前的原子集合体A进行模型化而获得的第1原子集合体模型。所制作的第1原子集合体模型的相关数据通过控制部12存储在存储部13中。在步骤S110中，控制部12起到第1原子集合体模型制作单元12A的作用。
在步骤S120中，控制部12利用存储在存储部13中的第1原子集合体模型的相关数据、第1模拟条件的相关数据13B、第1模拟程序13C等，通过针对第1原子集合体模型的计算机模拟，获取第1～第i的状态F1～Fi (i为2以上的整数。)的各状态下的原子集合体A的坐标。在步骤S120中，控制部12起到第1坐标获取单元12B的作用。
当采用分子动力学法作为计算机模拟时，控制部12所实行的步骤S120的一实施方式见图3。在该实施方式中，在步骤S110之后，控制部12对第1原子集合体模型实行能量最小化计算(步骤S121)，然后将第1原子集合体模型进行平衡化(步骤S122)，通过针对平衡化后的第1原子集合体模型的计算机模拟，获取平衡化后的状态F1～Fi的各状态下的原子集合体A的坐标。在步骤S122中，控制部12例如在第1原子集合体模型中的某物理量达到了阈值的阶段、或者从计算机模拟开始起经过了一定时间的阶段判断第1原子集合体模型达到了平衡化。
在步骤S130中，控制部12根据通过步骤S120获取的原子集合体A的坐标，获取通过片段B与状态F1～Fi的各状态下的原子集合体A结合而产生的原子集合体AB的坐标。在步骤S130中，控制部12起到第2坐标获取单元12C的作用。
控制部12所实行的步骤S130的一实施方式见图4。在该实施方式中，在步骤S120之后，控制部12实行以下的步骤：
步骤S131a，关于应该与状态Fx下的原子集合体A结合的片段B，规定选自构成片段B的原子的原子b相对于结构a的相对位置条件(例如相对于构成结构a的1个或2个以上的原子的距离、角度等)；
步骤S132a，根据所规定的相对位置条件，获取原子b的坐标；
步骤S133a，根据所获取的原子b的坐标和状态Fx下的原子集合体A的坐标RA(Fx)，获取通过原子b与状态Fx下的原子集合体A结合而产生的原子集合体Ab的坐标RAb(Fx)；以及
步骤S134a，以原子b的坐标为基准，利用计算机模拟、公知的构象产生方法等，使产生构成片段B的原子b以外的原子的坐标，获取通过片段B与状态Fx下的原子集合体A结合而产生的原子集合体AB的坐标RAB(Fx)。然后，在步骤S135a中，控制部12判定是否以x=1，2，···，i实行了步骤S131a～S134a，当为“NO”时重复步骤S131a～S134a。即，控制部12在步骤S135a中以x=1、2、···、i反复实行步骤S131a～S134a，直至判定为“YES”。由此，可以获取通过相对于结构a的相对位置已限定的片段B与状态F1～Fi的各状态下的原子集合体A结合而产生的原子集合体AB的坐标RAB(F1)～RAB(Fi)。
在步骤S131a中，关于应该与状态Fx下的原子集合体A结合的片段B，控制部12可以规定彼此不同的第1、第2、···、第p的相对位置条件。p是从x的值中独立选择的2以上的整数。因此，在x=1、2、···、i中p可以是相同的值，在x=1、2、···、i中p也可以是不同的值。
在步骤S131a中，关于应该与状态Fx下的原子集合体A结合的片段B，当规定彼此不同的第1、第2、···、第p的相对位置条件时，关于第1～第p的各相对位置条件，控制部12重复步骤S132a～S134a。由此，可以获取通过相对于结构a的相对位置不同的片段B1、B2、···、Bp与状态Fx下的原子集合体A结合而产生的原子集合体AB的坐标RAB(Fx，B1)、RAB(Fx，B2)、···、RAB(Fx，Bp)。
控制部12所实行的步骤S130的优选实施方式见图5。在该实施方式中，在步骤S120之后，控制部12实行以下的步骤：
步骤S131b，制作将原子集合体C进行模型化而获得的第3原子集合体模型，所述原子集合体C由结构a和与该结构a结合的片段B构成或者包含结构a和与该结构a结合的片段B而构成；以及
步骤S132b，通过针对所制作的第3原子集合体模型的计算机模拟，获取第1～第k的状态H1～Hk (k为2以上的整数。)的各状态下的原子集合体C的坐标；
之后，以x=1～i反复实行下述步骤：
步骤S133b，选择由选自构成结构a的原子的1个或2个以上的原子构成的选择原子组(需要说明的是，不仅包括选择2个以上的原子的情形，还包括选择1个原子的情形，为方便起见，称作“选择原子组”。)，使选自状态H1～Hk的1种或2种以上的状态下的原子集合体C的选择原子组的坐标相对于状态Fx下的原子集合体A的选择原子组的坐标旋转和/或并进，从而在原子集合体A的选择原子组与原子集合体C的选择原子组之间制作对应的原子间距离的平方和达到最小的原子集合体C的坐标(使原子集合体C的选择原子组的坐标与原子集合体A的选择原子组的坐标拟合)，根据所制作的原子集合体C的坐标，使选自状态H1～Hk的1种或2种以上的状态下的原子集合体C与原子集合体A重叠；以及
步骤S134b，根据状态Fx下的原子集合体A的坐标、以及与该原子集合体A重叠的原子集合体C的片段B的1个或2个以上的坐标，获取通过片段B与状态Fx下的原子集合体A结合而产生的原子集合体AB的坐标。
然后，在步骤S135b中，控制部12判定是否以x=1、2、···、i实行了步骤S133b～S134b，当为“NO”时重复步骤S133b～S134b。即，在步骤S135b中，控制部12以x=1、2、···、i反复实行步骤S133b～S134b，直至判定为“YES”。由此，可以获取通过相对于结构a的相对位置已限定的片段B与状态F1～Fi的各状态下的原子集合体A结合而产生的原子集合体AB的坐标RAB(F1)～RAB(Fi)。
在步骤S133b中，控制部12可以从状态H1～Hk下的原子集合体C的坐标中选择片段B相对于结构a的相对位置彼此不同的p组的坐标。p是从x的值中独立选择的2以上的整数。因此，在x=1、2、···、i中p可以是相同的值，在x=1、2、···、i中p也可以是不同的值。
在步骤S133b中，从状态H1～Hk下的原子集合体C的坐标中选择片段B相对于结构a的相对位置彼此不同的p组的坐标时，关于p组的各个坐标，控制部12重复步骤S133b～S134b。由此，可以获取通过相对于结构a的相对位置不同的片段B1、B2、···、Bp与状态Fx下的原子集合体A结合而产生的原子集合体AB的坐标RAB(Fx，B1)、RAB(Fx，B2)、···、RAB(Fx，Bp)。
在步骤S140中，控制部12根据通过步骤S130获取的原子集合体AB的坐标，计算结构a和与该结构a结合的片段B的相互作用能φ，在步骤S150中，制作表示相互作用能φ的各级频率的频率分布。在步骤S140和S150中，控制部12起到第1相互作用能φ频率分布制作单元12D的作用。
以下，对步骤S160～S190进行说明。
在步骤S150之后，控制部12实行相互作用能φ的相关处理(步骤S160)和相互作用能ε的相关处理(步骤S170～S190)。此时，控制部12可以在实行步骤S160之后再实行步骤S170～S190，也可以在实行步骤S170～S190之后再实行步骤S160，还可以并行实行步骤S160和步骤S170～S190。在任一种情况下，控制部12均依次实施步骤S170～S190。
在步骤S160中，控制部12根据通过步骤S150制作的频率分布，计算相互作用能φ的各级的出现概率P0(φ)。在步骤S160中，控制部12起到第1相互作用能φ出现概率计算单元12E的作用。
在步骤S170中，控制部12根据通过步骤S130获取的原子集合体AB的坐标，计算通过从原子集合体AB中去除由结构a和与该结构a结合的片段B构成的结构aB而产生的原子集合体的一部分或全部与片段B的相互作用能ε，在步骤S180中，制作通过步骤S150制作的频率分布的相互作用能φ的各级中的、表示相互作用能ε的各级频率的频率分布。在步骤S170和S180中，控制部12起到第1相互作用能ε频率分布制作单元12F的作用。
在步骤S170中，控制部12可以根据步骤S130中获取的原子集合体AB的坐标，关于原子集合体AB的各坐标RAB(F1)～RAB(Fi)，计算通过从原子集合体AB中去除由结构a和与该结构a结合的片段B构成的结构aB而产生的原子集合体的一部分或全部与片段B的相互作用能ε，之后，可以根据步骤S130中获取的原子集合体AB的坐标和步骤S150中制作的频率分布，抽取相互作用能φ的各级中的、通过从原子集合体AB中去除该结构aB而产生的原子集合体的一部分或全部与片段B的相互作用能ε，制作通过步骤S150制作的频率分布的相互作用能φ的各级中的、表示相互作用能ε的各级频率的频率分布，也可以根据步骤S150中制作的频率分布，从通过步骤S130获取的原子集合体AB的坐标中抽取属于相互作用能φ的各级的原子集合体AB的坐标，之后，根据所抽取的原子集合体AB的坐标，计算相互作用能φ的各级中的、通过从原子集合体AB中去除由结构a和与该结构a结合的片段B构成的结构aB而产生的原子集合体的一部分或全部与片段B的相互作用能ε，制作通过步骤S150制作的频率分布的相互作用能φ的各级中的、表示相互作用能ε的各级频率的频率分布。
在步骤S190中，控制部12根据步骤S180中制作的频率分布，计算相互作用能φ的各级中的、相互作用能ε的各级的出现概率P0’(ε；φ)。在步骤S190中，控制部12起到第1相互作用能ε出现概率计算单元12G的作用。
以下，对步骤S210～S240进行说明。
控制部12依次实施步骤S210～S240。
在步骤S210中，控制部12利用存储在存储部13中的变化后的原子集合体AB的相关数据13D，制作将变化后的原子集合体AB进行模型化而获得第2原子集合体模型。所制作的第2原子集合体模型的相关数据通过控制部12存储在存储部13中。在步骤S210中，控制部12起到第2原子集合体模型制作单元12H的作用。
在步骤S220中，控制部12利用存储在存储部13中的第2原子集合体模型的相关数据、第2模拟条件的相关数据13E、第2模拟程序13F等，通过针对第2原子集合体模型的计算机模拟，获取第1～第j的状态G1～Gj (j为2以上的整数。)的各状态下的原子集合体AB的坐标。在步骤S220中，控制部12起到第3坐标获取单元12I的作用。
当采用分子动力学法作为计算机模拟时，控制部12所实行的步骤S220的一实施方式见图6。在该实施方式中，在步骤S210之后，控制部12对第2原子集合体模型实行能量最小化计算(步骤S221)，然后将第2原子集合体模型进行平衡化(步骤S222)，通过针对平衡化后的第2原子集合体模型的计算机模拟，获取平衡化后的状态G1～Gj的各状态下的原子集合体AB的坐标。在步骤S222中，控制部12例如在第2原子集合体模型中的某物理量达到了阈值的阶段、或者从开始计算机模拟起经过了一定时间的阶段判断为第2原子集合体模型达到了平衡化。
在步骤S230中，控制部12根据步骤S220中获取的原子集合体AB的坐标，计算结构a和与该结构a结合的片段B的相互作用能φ，在步骤S240中，制作表示相互作用能φ的各级频率的频率分布。在步骤S230和S240中，控制部12起到第2相互作用能φ频率分布制作单元12J的作用。
以下，对步骤S250～S280进行说明。
在步骤S240之后，控制部12实行相互作用能φ的相关处理(步骤S250)和相互作用能ε的相关处理(步骤S260～S280)。此时，控制部12可以在实行步骤S250之后再实行步骤S260～S280，也可以在实行步骤S260～S280之后再实行步骤S250，还可以并行实行步骤S250和步骤S260～S280。在任一种情况下，控制部12均依次实施步骤S260～S280。
在步骤S250中，控制部12根据步骤S240中制作的频率分布，计算相互作用能φ的各级的出现概率P(φ)。在步骤S250中，控制部12起到第2相互作用能φ出现概率计算单元12K的作用。
在步骤S260中，控制部12根据步骤S220中获取的原子集合体AB的坐标，计算通过从原子集合体AB中去除由结构a和与该结构a结合的片段B构成的结构aB而产生的原子集合体的一部分或全部与片段B的相互作用能ε，在步骤S270中，制作表示相互作用能ε的各级频率的频率分布。在步骤S260和S270中，控制部12起到第2相互作用能ε频率分布制作单元12L的作用。
在步骤S260中，控制部12可以将相互作用能ε与步骤S240中制作的频率分布中的相互作用能φ的各级无关联地进行计算，也可以关联着进行计算。与相互作用能φ的各级关联着计算相互作用能ε时，相互作用能ε与步骤S240中制作的频率分布(例如直方图)中所含的相互作用能φ的所有级相关联进行计算，表示相互作用能ε的各级频率的频率分布与相互作用能φ的各级中的相互作用能ε相关联进行制作。
当控制部12将相互作用能ε与相互作用能φ的各级无关联地进行计算时，控制部12根据步骤S220中获取的原子集合体AB的坐标，关于原子集合体AB的各坐标RAB(G1)～RAB(Gj)，计算通过从原子集合体AB中去除由结构a和与该结构a结合的片段B构成的结构aB而产生的原子集合体的一部分或全部与片段B的相互作用能ε，制作表示相互作用能ε的各级频率的频率分布。
当控制部12将相互作用能ε与相互作用能φ的各级关联着进行计算时，控制部12根据步骤S220中获取的原子集合体AB的坐标，计算通过从原子集合体AB中去除由结构a和与该结构a结合的片段B构成的结构aB而产生的原子集合体的一部分或全部与片段B的相互作用能ε，制作通过步骤S240制作的频率分布的相互作用能φ的各级中的、表示相互作用能ε的各级频率的频率分布。
在一实施方式中，控制部12根据步骤S220中获取的原子集合体AB的坐标，关于原子集合体AB的各坐标RAB(G1)～RAB(Gj)，计算通过从原子集合体AB中去除由结构a和与该结构a结合的片段B构成的结构aB而产生的原子集合体的一部分或全部与片段B的相互作用能ε，之后，根据步骤S240中制作的频率分布，抽取相互作用能φ的各级中的、通过从原子集合体AB中去除由结构a和与该结构a结合的片段B构成的结构aB而产生的原子集合体的一部分或全部与片段B的相互作用能ε，制作通过步骤S240制作的频率分布的相互作用能φ的各级中的、表示相互作用能ε的各级频率的频率分布。
在又一实施方式中，控制部12根据步骤S240中制作的频率分布，从通过步骤S220获取的原子集合体AB的坐标中抽取属于相互作用能φ的各级的原子集合体AB的坐标，之后，根据所抽取的原子集合体AB的坐标，计算相互作用能φ的各级中的、通过从原子集合体AB中去除由结构a和与该结构a结合的片段B构成的结构aB而产生的原子集合体的一部分或全部与片段B的相互作用能ε，制作通过步骤S240制作的频率分布的相互作用能φ的各级中的、表示相互作用能ε的各级频率的频率分布。
当控制部12将相互作用能ε与相互作用能φ的各级无关联地进行计算时，在步骤S260和S270中，控制部12根据步骤S220中获取的原子集合体AB的坐标，关于原子集合体AB的各坐标RAB(G1)～RAB(Gj)，计算通过从原子集合体AB中去除由结构a和与该结构a结合的片段B构成的结构aB而产生的原子集合体的一部分或全部与片段B的相互作用能ε，制作表示相互作用能ε的各级频率的频率分布。
当控制部12将相互作用能ε与相互作用能φ的各级关联着进行计算时，在步骤S260和S270中，控制部12可以根据步骤S220中获取的原子集合体AB的坐标，关于原子集合体AB的各坐标RAB(G1)～RAB(Gj)，计算通过从原子集合体AB中去除由结构a和与该结构a结合的片段B构成的结构aB而产生的原子集合体的一部分或全部与片段B的相互作用能ε，之后，可以根据步骤S240中制作的频率分布，抽取相互作用能φ的各级中的、通过从原子集合体AB中去除结构aB而产生的原子集合体的一部分或全部与片段B的相互作用能ε，制作通过步骤S240制作的频率分布的相互作用能φ的各级中的、表示相互作用能ε的各级频率的频率分布，也可以根据步骤S240中制作的频率分布，从通过步骤S220获取的原子集合体AB的坐标中抽取属于相互作用能φ的各级的原子集合体AB的坐标，之后，根据所抽取的原子集合体AB的坐标，计算相互作用能φ的各级中的、从原子集合体AB中去除结构aB而获得的原子集合体的一部分或全部与片段B的相互作用能ε，制作通过步骤S240制作的频率分布的相互作用能φ的各级中的、表示相互作用能ε的各级频率的频率分布。
在步骤S280中，控制部12根据步骤S270中制作的频率分布，计算相互作用能ε的各级的出现概率P’(ε)。在步骤S280中，控制部12起到第2相互作用能ε出现概率计算单元12M的作用。
在步骤S260和S270中，当控制部12将相互作用能ε与步骤S240中制作的频率分布中的相互作用能φ的各级无关联地进行计算、并制作相互作用能ε的各级频率分布时，在步骤S280中，控制部12根据步骤S270中制作的频率分布，而与相互作用能φ的各级无关联，计算相互作用能ε的各级的出现概率P’(ε)。
在步骤S260和S270中，控制部12将相互作用能ε与通过步骤S240制作的频率分布中的相互作用能φ的各级关联着进行计算、并制作相互作用能φ的各级中的、相互作用能ε的各级频率分布时，在步骤S280中，控制部12根据步骤S270中制作的频率分布，计算相互作用能φ的各级中的、相互作用能ε的各级的出现概率P’(ε；φ)。
以下，对步骤S300进行说明。
控制部12在实行变化前的原子集合体A的相关计算处理(步骤S110～S190)和变化后的原子集合体AB的相关计算处理(步骤S210～S280)后，实行自由能变化量∫Δν(φ)P(φ)dφ的相关计算处理(步骤S300)。
在步骤S300中，控制部12根据步骤S250中计算的P(φ)、步骤S190中计算的P0’(ε；φ)和步骤S280中计算的P’(ε)，计算相互作用能ε所引起的自由能变化量∫Δν(φ)P(φ)dφ。在步骤S300中，控制部12起到∫Δν(φ)P(φ)dφ计算单元12N的作用。
在一实施方式中，控制部12不计算Δν(φ)而根据步骤S250中计算的P(φ)、步骤S190中计算的P0’(ε；φ)和步骤S280中计算的P’(ε)，利用能量表示法计算∫Δν(φ)P(φ)dφ。这里，使用的P’(ε)是指与相互作用能φ的各级无关联地进行计算的相互作用能ε的各级的出现概率。
在又一实施方式中，控制部12根据步骤S190中计算的P0’(ε；φ)和步骤S280中计算的P’(ε；φ)，计算相互作用能φ的各级中的、相互作用能ε所引起的自由能变化量Δν(φ)，根据所计算的Δν(φ)和步骤S250中计算的P(φ)来计算∫Δν(φ)P(φ)dφ。此时，∫Δν(φ)P(φ)dφ计算单元12N优选利用能量表示法，根据步骤S190中计算的P0’(ε；φ)和步骤S280中计算的P’(ε；φ)，计算相互作用能φ的各级中的、相互作用能ε所引起的自由能变化量Δν(φ)。
以下，对步骤S400进行说明。
控制部12在实行自由能变化量∫Δν(φ)P(φ)dφ的相关计算处理(步骤S300)之后，实行ΔG的相关计算处理(步骤S400)。
在步骤S400中，控制部12根据步骤S160中计算的P0(φ)、步骤S250中计算的P(φ)、步骤S300中计算的∫Δν(φ)P(φ)dφ和数学式(1)来计算ΔG，
[数学式11]
[式中，R表示气体常数，T表示反应式(1)所表示的变化发生时的绝对温度。]。
本实施方式所涉及的计算装置1的功能可以通过使用记录有程序(即，使控制部12实行步骤S110～S190、步骤S210～S280、步骤S300和步骤S400的程序)的计算机可读取的记录介质、并将该程序安装在计算机中来实现，所述程序使控制部12起到第1原子集合体模型制作单元12A、第1坐标获取单元12B、第2坐标获取单元12C、第1相互作用能φ频率分布制作单元12D、第1相互作用能φ出现概率计算单元12E、第1相互作用能ε频率分布制作单元12F、第1相互作用能ε出现概率计算单元12G、第2原子集合体模型制作单元12H、第3坐标获取单元12I、第2相互作用能φ频率分布制作单元12J、第2相互作用能φ出现概率计算单元12K、第2相互作用能ε频率分布制作单元12L、第2相互作用能ε出现概率计算单元12M、∫Δν(φ)P(φ)dφ计算单元12N和ΔG计算单元12O的作用。作为记录有程序的计算机可读取的记录介质，例如可以列举：ROM、软盘(注册商标)、硬盘、光盘、光磁盘、CD-ROM、磁带、非挥发性存储卡等。
需要说明的是，本发明不仅包含通过计算机执行程序来实现本实施方式所涉及的计算装置1的功能的实施方式，还包含通过程序与在计算机中运行的OS(操作系统)、其它应用软件等共同来实现本实施方式所涉及的计算装置1的功能的实施方式。另外，本发明还包含如下的实施方式：将程序存储在计算机的功能扩展板或与计算机连接的功能扩展单元所具备的存储器中，之后根据该程序的指示，该功能扩展板或该功能扩展单元所具备的CPU等实行实际处理的一部分或全部，从而实现本实施方式所涉及的计算装置1的功能。
实施例
以下，给出实施例，以阐述本发明的具体实施方式，但本发明并不受限于此。
在本实施例中，计算了在反应式(2)和(3)中以ΔG2-ΔG1定义的ΔΔG。
[化学式1]
在反应式(2)和(3)中，R表示蛋白，包围R的□表示在蛋白R的周围存在水分子，L1表示与蛋白R结合的配体，包围L1的□表示在配体L1的周围存在水分子，L1R表示蛋白R和与该蛋白R结合的配体L1的复合体，包围L1R的□表示在复合体L1R的周围存在水分子，L2是指与蛋白R结合的配体，表示不同于配体L1的配体，包围L2的□表示在配体L2的周围存在水分子，L2R表示蛋白R和与该蛋白R结合的配体L2的复合体，包围L2R的□表示在复合体L2R的周围存在水分子。需要说明的是，在反应式(2)的左边L1与R没有进行相互作用，而在反应式(3)的左边L2与R没有进行相互作用。
由反应式(4)可知：ΔΔG可以表示为ΔΔG=ΔG4-ΔG3。在反应式(4)的左边，蛋白R与配体L1和L2没有进行相互作用，所以ΔG3是指处于周围存在水分子的状态的配体L1与处于周围存在水分子的状态的配体L2的自由能差。
[化学式2]
根据反应式(5)，ΔG4可以表示为ΔG4=ΔG4a+ΔG4b。需要说明的是，在反应式(5)中，为方便起见，在L1R→L0R的变化中去除的片段和在L0R→L2R的变化中添加的片段均用“B”表示，但实际上在L1R→L0R的变化去除的片段与在L0R→L2R的变化中添加的片段不同。关于在L1R→L0R的变化中去除的片段和在L0R→L2R的变化中添加的片段的细节见后述。
[化学式3]
在反应式(5)中，L0R表示蛋白R和与该蛋白R结合的配体L0的复合体，包围L0R的□表示在复合体L0R的周围存在水分子。
根据反应式(6)，ΔG3可以表示为ΔG3=ΔG3a+ΔG3b。需要说明的是，在反应式(6)中，为方便起见，在L1→L0的变化中去除的片段和在L0→L2的变化中添加的片段均用“B”表示，但实际上在L1→L0的变化中去除的片段与在L0→L2的变化中添加的片段不同。关于在L1→L0的变化中去除的片段和在L0→L2的变化中添加的片段的细节见后述。
[化学式4]
在反应式(6)中，L0是指与蛋白R结合的配体，表示不同于配体L1和配体L2的配体，包围L0的□表示在配体L0的周围存在水分子。需要说明的是，在反应式(6)中，R与L1、L2、L0中的任一种配体均不发生相互作用。因此，反应式(6)中的ΔG3是指处于周围存在水分子的状态的配体L1与处于周围存在水分子的状态的配体L2的自由能差。
[实施例1]
(1) ΔG4a的计算
关于反应式(7)所表示的变化，以复合体L0R和构成存在于该复合体L0R周围的水分子的变化前的原子集合体作为“原子集合体A”、以复合体L1R和构成存在于该复合体L1R周围的水分子的变化后的原子集合体作为“原子集合体AB”，根据数学式(1)计算了变化前的原子集合体A的自由能和片段B的自由能的和、与变化后的原子集合体AB的自由能之差ΔG4a’。
[化学式5]
需要说明的是，由于ΔG4a是与反应式(8)所表示的变化有关的自由能变化量，因此在根据数学式(1)计算的ΔG4a’上乘以负号(minus)而得到的值相当于ΔG4a。
[化学式6]
作为蛋白R，选择了作为丝氨酸蛋白酶之一的胰蛋白酶(Trypsin)。
作为配体L1，选择了作为胰蛋白酶抑制剂的苯甲脒。配体L1的化学结构如下。
[化学式7]
构成配体L1的各原子的部分电荷如下。关于部分电荷，将小数点以下第4位四舍五入，记载到小数点以下第3位。以后所记载的部分电荷亦同。需要说明的是，这些部分电荷是通过后述的MATCH程序制作的。
[化学式8]
作为配体L0，选择了从配体L1中去除片段B而获得的分子，所述片段B由位于脒基对位的苯环的碳原子上结合的氢原子和具有该氢原子所结合的苯环的碳原子的部分电荷(-0.115)的点电荷构成。配体L0的化学结构如下。
[化学式9]
构成配体L0的各原子的部分电荷如下。
[化学式10]
作为片段B，选择了氢原子和点电荷(部分电荷；-0.115)。氢原子是指在配体L1中与苯环结合的氢原子，点电荷在配体L1中具有结合有氢原子的碳原子所具有的部分电荷(-0.115)。通过从配体L1中去除片段B，获得在构成苯环的碳原子中片段B的氢原子所结合的碳原子的部分电荷为零的分子作为配体L0。片段B的结构如下。这里，具有部分电荷(-0.115)的点电荷与氢原子并未结合，即不存在由键长、键角、扭转角等所引起的相互作用。因此，以下虚拟原子与氢原子之间用虚线进行图示。另外，片段B的点电荷与点电荷以外的片段B的原子的静电相互作用，假定在该原子与点电荷之间具有共价键时，在该原子间包含2个以上的其它原子的情况下产生。在以下的其它实施例中也同样记载。
[化学式11]
从配体L0和配体L1的结构中选择构成配体L0的所有原子，作为“结构a”。因此，配体L0由结构a构成，而配体L1由结构a和与该结构a结合的片段B构成。
在本实施例中，利用图1所示的计算装置，实行图2所示的步骤S110～S190、S210～S280、S300和S400所对应的步骤S110A～S190A、S210A～S280A、S300A和S400A，从而算出了ΔG4a’。在步骤S120A中实行了图3所示的步骤S121～S123所对应的步骤S121A～S123A，在步骤S130A中实行了图5所示的步骤S131b～S134b所对应的步骤S131A～S134A，在步骤S220A中实行了图6所示的步骤S221～S223所对应的步骤S221A～S223A。
[变化前的原子集合体A所包含的复合体L0R和变化后的原子集合体AB所包含的复合体L1R的相关数据的输入]
作为实行步骤S110A和S210A之前的准备步骤，输入部11将变化前的原子集合体A所包含的复合体L0R和变化后的原子集合体AB所包含的复合体L1R的相关数据(构成复合体L1R的各原子的坐标、种类、质量、部分电荷、原子间的结合信息等，构成复合对L0R的各原子的坐标、种类、质量、部分电荷、原子间的结合信息等)输入到控制部12中。这里，原子的种类是指元素周期表中的元素的种类和原子类型，各原子所具有的原子类型或各原子间的原子类型对是指与势能参数的值一对一地对应赋予的类型。所输入的数据通过控制部12存储在存储部13中。
变化后的原子集合体AB所包含的复合体L1R的相关数据是利用立体结构数据文件(3PTB_AB)和可视化软件VMD(Visual Molecular Dynamics：http：//www.ks.uiuc.edu/Research/vmd/)中安装的“Automatic PSF Builder”功能得到的，所述立体结构数据文件(3PTB_AB)如下制作：利用由CCG(Chemical Computing Group)公司提供的综合计算化学系统MOE(Molecular Operating Environment)，相对于从作为立体结构数据库的蛋白数据库(Protein Data Bank，PDB；http：//www.rcsb.org/pdb/home/home.do)中下载的立体结构数据文件(PDB code：3PTB)，删除该立体结构数据文件所包含的归属于水分子的氧原子，之后利用MOE中安装的“Protonate 3D”功能追加氢原子和该氢原子的坐标而制作立体结构数据文件(3PTB_AB)。需要说明的是，立体结构数据(PDB code：3PTB)包含与钙有关的坐标数据，该钙作为离子、还作为蛋白的一部分来处理。另外，变化前的原子集合体A所包含的复合体L0R的相关数据是根据立体结构数据文件(3PTB_A)，利用VMD的“Automatic PSFBuilder”功能制作的，所述立体结构数据文件(3PTB_A)是利用MOE的“Builder”功能，从立体结构数据文件(3PTB_AB)中去除作为构成片段B的虚拟原子的点电荷以外的原子而制作的。
需要说明的是，在立体结构数据文件(3PTB_AB)中包含复合体L1R的相关坐标数据，利用MOE从立体结构数据文件(3PTB_AB)中选择配体L1的相关坐标数据，制作配体L1的相关立体结构数据文件(3PTB_L1)，利用VMD的“Automatic PSF Builder”功能，可以获得配体L1的相关数据。同样，由立体结构数据文件(3PTB_A)制作配体L0的相关立体结构数据文件(3PTB_L0)，并且，关于部分电荷，将与构成去除了虚拟原子的片段B的原子(即氢原子)结合的苯环上的碳原子的部分电荷设为零，从而得到了配体L0的相关数据。
[模拟条件的相关数据的输入]
在后述的实行步骤S110A和S210A之前的准备步骤中，输入部11将能量最小化计算的相关模拟条件(对各原子施加的力的计算次数、势能参数的种类和值、边界条件、范德华势能和库仑势能计算的相关转换函数和截断半径、库仑势能计算中的长距离相互作用、1-4相互作用的相关条件等)和分子动力学模拟的相关模拟条件(模拟时间、温度条件、压力条件、势能参数的种类和值、产生的系综的种类、边界条件、运动方程式的数值解法、数值积分的时间间隔(時間刻み)、范德华势能和库仑势能计算的相关转换函数和截断半径、库仑势能计算中的长距离相互作用的相关条件、1-4相互作用、抽点打印数目等的输出条件等)分别输入到控制部12中。所输入的数据通过控制部12存储在存储部13中。
构成蛋白R的各原子的势能参数的种类使用了CHARMm22，构成配体L0的各原子的势能参数的种类和构成配体L1的各原子的势能参数的种类使用了CGenFF(Charmm GeneralForceField)，构成水分子的各原子的势能参数的种类使用了TIP3P。需要说明的是，配体L0的各原子的势能参数和配体L1的各原子的势能参数利用MATCH程序(http：//brooks.chem.lF.umich.edu/index.phppage=match＆subdir=articles/resources/soFtware)进行了确定。
[步骤S110A：第1原子集合体模型的制作]
在步骤S110A中，控制部12利用可视化软件VMD读取存储在存储部13中的立体结构数据文件(3PTB_A)，然后，利用该VMD中安装的“solvate”功能，制作使在蛋白R和与该蛋白R结合的配体L0的复合体L0R的周围产生了9492个水分子的立方体型基本单元，制作将变化前的原子集合体A进行模型化而获得的第1原子集合体模型。所制作的第1原子集合体模型的相关数据(构成原子集合体A的各原子的坐标、种类、质量、部分电荷、原子间的结合信息等)通过控制部12存储在存储部13中。这里，原子的种类是指元素周期表中的元素的种类和指定用于计算机模拟的势能参数的原子类型。
[步骤S120A：通过计算机模拟获取原子集合体A的坐标]
在步骤S120A中，控制部12实行以下的步骤S121A～S123A。
[步骤S121A：针对第1原子集合体模型的能量最小化计算]
在步骤S121A中，对于第1原子集合体模型，控制部12将能量最小化计算的相关模拟条件(对各原子施加的力的计算次数、势能参数的种类和值、边界条件、范德华势能和库仑势能计算的相关转换函数和截断半径、库仑势能计算中的长距离相互作用、1-4相互作用的相关条件等)读取到作为分子动力学模拟程序的NAMD程序(http：//www.ks.uiuc.edu/Research/namd/)中，实施利用了NAMD的“minimize”功能的能量最小化计算。能量最小化计算后的第1原子集合体模型的相关数据(构成原子集合体A的各原子的坐标、种类、质量、部分电荷、原子类型、原子间的结合信息等)通过控制部12存储在存储部13中。这里，原子的种类是指元素周期表中的元素的种类和指定计算机模拟中使用的势能参数的原子类型。
能量最小化计算的条件如下。
采用周期边界条件作为边界条件，关于计算次数，在NAMD程序的输入文件中，将“minimize(最小化)”设为10000。
关于转换函数和截断半径，在NAMD程序的输入文件中，将用于计算范德华势能和库仑势能的转换函数的相关关键词“switching”设为10埃，将截断半径的相关关键词“cutoff”设为12埃。
关于库仑势能计算的长距离相互作用，采用了将“PMEGridSpacing”设为1.0埃和将“PMEInterpOrder”设为4的颗粒网状埃瓦尔德(PME：Particle Mesh Ewald)法。
关于1-4相互作用，将“exclude”设为scaled1-4、将“1-4scaling”设为1。
[步骤S122A：第1原子集合体模型的平衡化]
在步骤S122A中，控制部12将存储在存储部13中的能量最小化计算后的第1原子集合体模型的相关数据(构成原子集合体A的各原子的坐标、种类、质量、部分电荷、原子间的结合信息等)、模拟条件的相关数据(模拟时间、温度条件、压力条件、势能参数的种类和值、产生的系综的种类、边界条件、运动方程式的数值解法、数值积分的时间间隔、关于计算范德华势能和库仑势能的转换函数和截断半径、库仑势能计算中的长距离相互作用、1-4相互作用的相关条件、抽点打印数目等的输出条件等)读取到NAMD程序中，对能量最小化计算后的第1原子集合体模型实行分子动力学模拟，使第1原子集合体模型进行了平衡化。平衡化后的第1原子集合体模型的相关数据(构成原子集合体A的各原子的坐标、种类、质量、部分电荷、原子间的结合信息等)通过控制部12存储在存储部13中。这里，原子的种类是指元素周期表中的元素的种类和指定用于计算机模拟的势能参数的原子类型。
需要说明的是，分子动力学模拟的条件如下进行了设定。
采用周期边界条件作为边界条件。关于压力条件，在NAMD程序的输入文件中，将用于控制压力的关键词“langevinPiston”设定为on、并将“langevinPistonTarget”设定为1.01325，从而将压力控制在1atm。关于温度条件，在NAMD程序的输入文件中，将用于控制温度的关键词“langevin”设定为on、并将“langevinTemp”以每50开尔文进行设定直至达到50开尔文～300开尔文，从而使之阶段性升温。模拟时间在50开尔文～250开尔文下分别设为100皮可秒，在300开尔文下设为1000皮可秒。作为运动方程式的数值解法，采用shake法，将时间间隔设为2飞秒。作为输出条件，使每200飞秒输出第1原子集合体模型的坐标。
关于转换函数和截断半径，在NAMD程序的输入文件中，将用于计算范德华势能和库仑势能的转换函数的相关关键词“switching”设定为10埃，将截断半径的相关关键词“cutoff”设定为12埃。
关于库仑势能的计算的长距离相互作用，采用了将“PMEGridSpacing”设定为1.0埃和将“PMEInterpOrder”设定为4的颗粒网状埃瓦尔德(PME：Particle Mesh Ewald)法。
关于1-4相互作用，将“exclude”设为scaled1-4、将“1-4scaling”设为1。
关于系综，通过实行安装在NAMD程序中上述温度控制和压力控制，构成了NPT系综。
[步骤S123A：通过计算机模拟获取原子集合体A的坐标]
在步骤S123A中，控制部12将存储在存储部13中的平衡化后的第1原子集合体模型的相关数据(构成原子集合体A的各原子的坐标、种类、质量、部分电荷、原子间的结合信息等)、模拟条件的相关数据(模拟时间、温度条件、压力条件、势能参数的种类和值、产生的系综的种类、边界条件、运动方程式的数值解法、数值积分的时间间隔、关于计算范德华势能和库仑势能的转换函数和截断半径、库仑势能计算中的长距离相互作用的相关条件、1-4相互作用抽点打印数目等的输出条件等)读取到NAMD程序中，对平衡化后的第1原子集合体模型实行了分子动力学模拟。这里，原子的种类是指元素周期表中的元素的种类和指定用于计算机模拟的势能参数的原子类型。
需要说明的是，分子动力学模拟的条件如下进行了设定。
采用周期边界条件作为边界条件。关于压力条件，在NAMD程序的输入文件中，将用于控制压力的关键词“langevinPiston”设定为on、将“langevinPistonTarget”设定为1.01325，从而将压力控制在1atm。关于温度条件，在NAMD程序的输入文件中，将用于控制温度的关键词“langevin”设定为on、将“langevinTemp”设定为300开尔文。模拟时间设为5000皮可秒，采用shake法作为运动方程式的数值解法，时间间隔设为2飞秒。关于输出条件，使每50飞秒输出第1原子集合体模型的坐标。
关于转换函数和截断半径，在NAMD程序的输入文件中，将用于计算范德华势能和库仑势能的转换函数的关键词“switching”设为10埃、将关于截断半径的关键词“cutoff”设为12埃。
关于库仑势能的计算的长距离相互作用，采用了将“PMEGridSpacing”设定为1.0埃和将“PMEInterpOrder”设定为4的颗粒网状埃瓦尔德(PME：Particle Mesh Ewald)法。
关于1-4相互作用，将“exclude”设为scaled1-4、将“1-4scaling”设为1。
关于系综，通过实行安装在NAMD程序中的上述温度控制和压力控制，构成了NPT系综。
控制部12从开始分子动力学模拟起每50飞秒获取了即50飞秒后(时间T1)、100飞秒后(时间T2)、···、5000皮可秒后(时间T100000)的原子集合体A的坐标。由此，获取了时间Ti(i=1，2，···，100000)时的状态Fi下的原子集合体A的坐标RA(Fi)。所获取的原子集合体A的坐标为RA(F1)、RA(F2)、···、RA(F100000)的总计100000组的坐标。
控制部12将原子集合体A的坐标RA(F1)～RA(F100000)沿时间轴排列存储在存储部13中。即，控制部12将时间Ti(i=1，2，···，100000)时的状态Fi下的原子集合体A的坐标RA(Fi)与时间Ti对应记入存储在存储部13中。
[步骤S130A：获取原子集合体AB的坐标]
在步骤S130A中，控制部12实行了步骤S131A～S134A。
[步骤S131A：第3原子集合体模型的制作]
选择了由结构a和与该结构a结合的片段B构成的配体L1作为原子集合体C。在原子集合体C中，配体L1的周围为真空。
在步骤S131A中，控制部12将存储在存储部13中的第3原子集合体模型的相关数据(构成配体L1的各原子的坐标、种类、质量、部分电荷、原子间的结合信息等)、关于能量最小化计算的模拟条件(对各原子施加的力的计算次数、势能参数的种类和值、边界条件、关于计算范德华势能和库仑势能的转换函数和截断半径、库仑势能计算中的长距离相互作用、1-4相互作用的相关条件等)的相关数据读取到NAMD程序中，对第3原子集合体模型实施了能量最小化计算。能量最小化计算后的第3原子集合体模型的相关数据(构成配体L1的各原子的坐标、种类、质量、部分电荷、原子间的结合信息等)通过控制部12存储在存储部13中。这里，原子的种类是指元素周期表中的元素的种类和指定用于计算机模拟的势能参数的原子类型。
需要说明的是，能量最小化计算的条件如下。
关于边界条件，在步骤S131A中，由于是真空中的配体L1的能量最小化计算，所以没有特别指定边界条件和库仑势能计算的长距离相互作用的相关条件。关于对各原子施加的力的计算次数，在NAMD程序的输入文件中，将“minimize”设为1000。
关于转换函数和截断半径，在NAMD程序的输入文件中，将用于计算范德华势能和库仑势能的转换函数的相关关键词“switching”设定为10埃、将关于截断半径的关键词“cutoff”设定为12埃。
关于1-4相互作用，将“exclude”设为scaled1-4、将“1-4scaling”设为1。
[步骤S132A：通过计算机模拟获取原子集合体C的坐标]
在步骤S132A中，控制部12对能量最小化计算后的第3原子集合体模型实行了分子动力学模拟。
分子动力学模拟条件如下进行了设定。
关于温度条件，在NAMD程序的输入文件中，将用于控制温度的关键词“langevin”设定为on，将“langevinTemp”设定为300开尔文。模拟时间设为10纳秒，采用shake法作为运动方程式的数值解法，时间间隔设为2飞秒。关于输出条件，每隔100飞秒输出第3原子集合体模型的坐标，使生成100000组的原子集合体C的坐标。
关于转换函数和截断半径，在NAMD程序的输入文件中，将用于计算范德华势能和库仑势能的转换函数的相关关键词“switching”设定为10埃、将关于截断半径的关键词“cutoff”设定为12埃。
关于1-4相互作用，将“exclude”设为scaled1-4、将“1-4scaling”设为1。
控制部12从开始分子动力学模拟起每100飞秒获取了即100飞秒后(时间T1)、200飞秒后(时间T2)、····、10纳秒后(时间T100000)的原子集合体C的坐标。由此，获取了时间Tk(k=1、2、···、100000)所对应的状态Hk下的原子集合体C的坐标RC(Hk)。所获取的原子集合体C的坐标为RC(H1)、RC(H2)、···、RC(H100000)的总计100000组的坐标。
控制部12将原子集合体C的坐标RC(H1)～RC(H100000)沿时间轴排列存储在存储部13中。即，控制部12将时间Tk(k=1、2、···、100000)所对应的状态Hk下的原子集合体C的坐标RC(Hk)与时间Tk对应记入存储在存储部13中。
[步骤S133A：状态Fx下的原子集合体C的坐标与原子集合体A的坐标的拟合]
在步骤S133A中，控制部12构成了原子集合体A的坐标RA(F1)与原子集合体C的坐标RC(H1)的对Q(F1，H1)、原子集合体A的坐标RA(F2)与原子集合体C的坐标RC(H2)的对Q(F2，H2)、····、原子集合体A的坐标RA(F100000)与原子集合体C的坐标RC(H100000)的对Q(F100000，H100000)的总计100000组的对。
在各对中，选择构成结构a(配体L0)的所有原子，作为“选择原子组”。
关于对Q(F1，H1)～Q(F100000，H100000)的各对，控制部12在保持原子集合体C的坐标的各原子间的相对坐标(内部坐标)的同时，利用数学式(7)，使原子集合体C中的构成选择原子组的原子与原子集合体A中的构成选择原子组的原子拟合，从而使原子集合体C与原子集合体A重叠。此时，在数学式(7)中针对各原子而设定的Cn如下。
[化学式12]
[步骤S134A：原子集合体AB的坐标的获取]
在步骤S134A中，关于对Q(F1，H1)～Q(F100000，H100000)的各对，控制部12从与原子集合体A重叠的原子集合体C的坐标中获取原子集合体C所具有的片段B的坐标(其中，作为虚拟原子的点电荷的坐标除外)，将其追加到原子集合体A的坐标中，获取原子集合体AB的坐标RAB(Fi)。由此，获取了原子集合体AB的坐标RAB(F1)～RAB(F100000)。
[步骤S140A：相互作用能φ的计算]
在步骤S140A中，关于原子集合体AB的各坐标RAB(F1)～RAB(F100000)，控制部12计算了结构a和与该结构a结合的片段B的相互作用能φ。此时，控制部12实行了以下的步骤。
控制部12以i=1、2、···、100000反复实行从原子集合体AB的坐标RAB(Fi)中抽取构成由结构a和与该结构a结合的片段B构成的结构aB(配体L1)的各原子的坐标RL1(Fi)的处理，获取了从原子集合体AB的坐标RAB(F1)～RAB(F100000)中抽取的、构成结构aB(配体L1)的各原子的坐标RL1(F1)～RL1(F100000)。
控制部12将构成结构aB(配体L1)的各原子的坐标、种类、质量、部分电荷、势能参数的种类和值、原子间的结合信息等读取到NAMD中，以i=1、2、···、100000反复实行计算坐标RL1(Fi)中的结构aB(配体L1)的能EL1(Fi)的处理，关于各坐标RL1(F1)～RL1(F100000)，计算了结构aB(配体L1)的能EL1(F1)～EL1(F100000)。这里，EL1(Fi)是结构aB(配体L1)的原子间的相互作用所引起的能，是配体L1的内部能量。
控制部12以i=1、2、···、100000反复实行从原子集合体AB的坐标RAB(Fi)中抽取构成结构a(配体L0)的各原子的坐标Ra(Fi)的处理，获取了从原子集合体AB的坐标RAB(F1)～RAB(F100000)中抽取的、构成结构a的各原子的坐标Ra(F1)～Ra(F100000)。
控制部12使与构成结构a的各原子有关的坐标、种类、质量、部分电荷、势能参数的种类和值、原子间的结合信息等读取到NAMD中，以i=1、2、···、100000反复实行计算坐标Ra(Fi)中的结构a的能Ea(Fi)的处理，关于各坐标Ra(F1)～Ra(F100000)，计算了结构a的能Ea(F1)～Ea(F100000)。这里，Ea(Fi)是结构a的原子间的相互作用所引起的能量，是结构a的内部能量。
控制部12以i=1、2、···、100000反复实行从原子集合体AB的坐标RAB(Fi)中获取构成去除了虚拟原子的片段B的原子(在结构aB为配体L1的本实施例中是指氢原子)的坐标Rd(Fi)的处理，获取了从原子集合体AB的坐标RAB(F1)～RAB(F100000)中抽取的、构成去除了虚拟原子的片段B的原子的坐标Rd(F1)～Rd(F100000)。
控制部12使与构成去除了虚拟原子的片段B的原子相关的坐标、种类、质量、部分电荷、势能参数的种类和值、原子间的结合信息等读取到NAMD中，关于坐标Rd(Fi)，以i=1、2、···、100000反复实行计算构成去除了虚拟原子的片段B的原子的能Ed(Fi)和库仑势能ECd(Fi)的处理，关于各坐标Rd(F1)～Rd(F100000)，获取了构成去除了虚拟原子的片段B的原子的能Ed(F1)～Ed(F100000)和库仑势能ECd(F1)～ECd(F100000)。这里，Ed(Fi)是构成去除了虚拟原子的片段B的原子间的相互作用所引起的能，ECd(Fi)是构成去除了虚拟原子的片段B的原子间的静电相互作用所引起的能。需要说明的是，当结构aB为配体L1时，构成去除了虚拟原子的片段B的原子仅为1个氢原子，因此所有的值均为零。
关于坐标Rd(Fi)，控制部12以i=1、2、···、100000反复实行从构成去除了虚拟原子的片段B的原子的能Ed(Fi)中减去库仑势能ECd(Fi)以计算能Ee(Fi)的处理，计算了能Ee(F1)～Ee(F100000)。
控制部12以i=1、2、···、100000反复实行由原子集合体AB的坐标RAB(Fi)获取构成结构f的各原子(即，构成包含虚拟原子的片段B的原子)的坐标Rf(Fi)的处理，所述结构f由构成去除了虚拟原子的片段B的原子和在构成结构aB(配体L1)的原子中具有构成片段B的点电荷的部分电荷的原子(在构成配体L1的苯环的碳原子中，与构成片段B的氢原子具有键的碳原子)构成，获取了从原子集合体AB的坐标RAB(F1)～RAB(F100000)中抽取的、构成结构f的各原子的坐标Rf(F1)～Rf(F100000)。
控制部12使构成结构f的各原子的坐标、种类、质量、部分电荷、势能参数的种类和值、原子间的结合信息等读取到NAMD中，关于坐标Rf(Fi)，以i=1、2、···、100000反复实行计算库仑势能ECf(Fi)的处理，获取了库仑势能ECf(F1)～ECf(F100000)。这里，ECf(Fi)是构成结构f的原子间的静电相互作用所引起的能。需要说明的是，在结构aB为配体L1的本实施例中，当假定在构成片段B的氢原子与点电荷之间具有共价键时，该原子间不包含其它原子，所以ECf(Fi)为零。
控制部12以i=1、2、···、100000反复实行补充能Ee(Fi)和库仑势能ECf(Fi)来计算能EB(Fi)的处理，计算了能EB(F1)～EB(F100000)。
控制部12根据下式，以i=1、2、···、100000反复实行了计算结构a与片段B的相互作用能φ的处理。
[步骤S150A：相互作用能φ的各级频率分布的制作]
在步骤S150A中，控制部12制作了直方图，其中横轴表示相互作用能φ的各级，纵轴表示相互作用能φ的各级频率。需要说明的是，在与相互作用能φ相关的直方图中，对等级间隔Δφ设定各级间隔Δφ，使各级间隔Δφ中的该Δφ与其频率H0(φ)之积达到恒定，相互作用能φ的分配函数、即等级间隔Δφ的个数设为100。
[步骤S160A：相互作用能φ的各级的出现概率P0(φ)的计算]
在步骤S160A中，控制部12将步骤S150A中制作的直方图标准化，计算了相互作用能φ的各级的出现概率P0(φ)。
[步骤S170A：相互作用能ε的计算]
输入部11将利用作为能量表示法程序组(http：//sourceForge.net/p/ermod/wiki/Home/)的一部分的gen_structure、gen_input程序制作的参数文件等计算所需的各种文件输入到控制部12中。所输入的文件通过控制部12存储在存储部13中。需要说明的是，本发明中使用的能量表示法程序组采用了0.2.3版。
需要说明的是，关于利用gen_input程序制作的、输入有构成结构aB(配体L1)的各原子的部分电荷和范德华的相关参数的SltInfo文件，将构成结构a的原子的范德华的相关参数变为零。另外，将构成包含虚拟原子的片段B的原子的部分电荷以外也变为零。
控制部12利用作为能量表示法程序组的一部分的ermod程序，关于原子集合体AB的坐标RAB(Fi)，以i=1、2、···、100000反复实行了计算片段B与蛋白R(本实施例中为1个)的相互作用能εp(i)、以及片段B与各水分子(在本实施例中为9492个的各水分子)的相互作用能εw(i)的处理。这里，各εw(i)是指具有片段B与各水分子的总计9492个相互作用能的值的集合。
[步骤S180A：相互作用能φ的各级中的、相互作用能ε的各级频率分布的制作]
控制部12从步骤S170A中计算的εp(1)～εp(100000)和εw(1)～εw(100000)中选择：利用步骤S150A中制作的直方图中的属于相互作用能φ的各级的坐标计算的εp和εw，制作了与所选择的εp和εw相关的直方图H0’(εp；φ)、H0’(εw；φ)。
[步骤S190A：相互作用能φ的各级中的、相互作用能ε的各级的出现概率P0’(ε；φ)的计算]
控制部12计算了直方图H0’(εp；φ)和H0’(εw；φ)中的各级的出现概率P0’(εp；φ)和P0’(εw；φ)。
这里，步骤S180A和S190A中的、直方图H0’(εp；φ)、H0’(εw；φ)和出现概率P0’(εp；φ)、P0’(εw；φ)通过作为能量表示法程序组的一部分的ermod程序来制作。需要说明的是，关于ermod程序的参数，在直方图H0’(εp；φ)和出现概率P0’(εp；φ)的制作中，在输入ermod程序参数的parameter_er文件中，根据缺省条件将ecfbin和ec0bin的设定变更为0.05、将engdiv的设定变更为10、将maxins的设定变更为1。另外，根据缺省条件将engdiv的设定变更为10、将maxins的设定变更为1，实施H0’(εw；φ)和P0’(εw；φ)的制作。
控制部12对步骤S150A中制作的直方图中的相互作用能φ的各级实施步骤S180A～S190A，关于相互作用能φ的各级，计算了相互作用能ε的出现概率P0’(εp；φ)和P0’(εw；φ)。
[步骤S210A：第2原子集合体模型的制作]
[步骤S220A：通过计算机模拟获取原子集合体AB的坐标]
在步骤S210A和S220A中，除了将原子集合体A变更为原子集合体AB以外，控制部12进行与步骤S110A和S120A相同的操作，实行步骤S210A和S220A，从开始步骤S223A的分子动力学模拟起，获取了每50飞秒、即50飞秒后(时间T1)、100飞秒后(时间T2)、····、5000皮可秒后(时间T100000)的原子集合体AB的坐标。由此，获取了时间Tj(j=1，2，···，100000)时的状态Gj下的原子集合体AB的坐标RAB(Gj)。所获取的原子集合体AB的坐标为RAB(G1)、RAB(G2)、···、RAB(G100000)的总计100000组的坐标。
控制部12将原子集合体AB的坐标RAB(G1)～RAB(G100000)沿时间轴排列存储在存储部13中。即，控制部12将时间Tj(j=1，2，···，100000)时的状态Gj下的原子集合体AB的坐标RAB(Gj)与时间Tj对应记入存储在存储部13中。
[步骤S230A：相互作用能φ的计算]
在步骤S230A中，控制部12进行与步骤S140A相同的操作，实行步骤S230A，关于原子集合体AB的各坐标RAB(G1)～RAB(G100000)，计算了结构a和与该结构a结合的片段B的相互作用能φ。
[步骤S240A：相互作用能φ的各级频率分布的制作]
在步骤S240A，控制部12制作了直方图，其中横轴表示相互作用能φ的各级，纵轴表示相互作用能φ的各级频率。各级间隔Δφ设为与步骤S150A中设定的各Δφ相同的Δφ，相互作用能φ的分配函数、即等级间隔Δφ的个数设为100。
[步骤S250A：相互作用能φ的各级的出现概率P(φ)的计算]
在步骤S250A中，控制部12将步骤S240A中制作的直方图进行标准化，计算了相互作用能φ的各级的出现概率P(φ)。
[步骤S260A：相互作用能ε的计算]
在步骤S260A中，控制部12利用作为能量表示法程序组的一部分的ermod程序，关于原子集合体AB的坐标RAB(Gj)，以i=1、2、···、100000反复实行了计算片段B与蛋白R(本实施例中为1个)的相互作用能εp(i)、以及片段B与各水分子(本实施例中为9492个的各水分子)的相互作用能εw(i)的处理。这里，各εw(i)是具有片段B与各水分子的总计9492个的相互作用能的值的集合。此时，进行与步骤S170A相同的操作，制作计算所需的各种文件，同时变更了SltInfo文件。
[步骤S270A：相互作用能ε的各级频率分布的制作]
在步骤S270A中，控制部12制作了与步骤S260A中计算的εp(1)～εp(100000)和εw(1)～εw(100000)相关的直方图H’(εp)、H’(εw)。
[步骤S280A：相互作用能ε的各级的出现概率P’(ε)的计算]
在步骤S280A中，控制部12计算了直方图H’(εp)和H’(εw)中的各级的出现概率P’(εp)和P’(εw)。
这里，步骤S270A和S280A中的直方图H’(εp)、H’(εw)和出现概率P’(εp)、P’(εw)通过作为能量表示法程序组的一部分的ermod程序来制作。关于ermod程序的参数，在直方图H’(εp)和P’(εp)的制作中，在输入ermod程序参数的parameter_er文件中，根据缺省条件将ecfbin和ec0bin的设定变更为0.05。另外，H’(εw)和P’(εw)的制作按照缺省条件来实施。
控制部12通过步骤S260A～S280A计算了相互作用能ε的各级的出现概率P’(εp)和P’(εw)。
[步骤S300A：相互作用能ε所引起的自由能变化量∫Δν(φ)P(φ)dφ的计算]
在步骤S300A中，控制部12通过以下的步骤，计算了相互作用能ε所引起的自由能变化量∫Δν(φ)P(φ)dφ。
控制部12利用步骤S250A中计算的相互作用能φ的各级的出现概率P(φ)、步骤S190A中计算的相互作用能φ的各级中的相互作用能εp的各级的出现概率P0’(εp：φ)、步骤S280A中计算的相互作用能εp的各级的出现概率P’(εp)、以及作为能量表示法程序组的一部分的slvfe程序，计算了∫Δνp(φ)P(φ)dφ。
控制部12利用步骤S250A中计算的相互作用能φ的各级的出现概率P(φ)、步骤S190A中计算的相互作用能φ的各级中的相互作用能εw的各级的出现概率P0’(εw：φ)、步骤S280A中计算的相互作用能εw的各级的出现概率P’(εw)、以及作为能量表示法程序组的一部分的slvfe程序，计算了∫Δνw(φ)P(φ)dφ。
控制部12利用作为slvfe程序的输出而得到的片段B的自身相互作用能Eself(Self-energy)和上述中计算的∫Δνp(φ)P(φ)dφ和∫Δνw(φ)P(φ)dφ，通过数学式(9)，计算了相互作用能ε所引起的自由能变化量∫Δν(φ)P(φ)dφ。需要说明的是，Eself是指：在能量表示法中，在考虑到周期边界条件而实施分子动力学模拟时，作为修正项而算出的值。
[步骤S400A：ΔG4a的计算]
在步骤S400A中，控制部12根据步骤S160A中计算的相互作用能φ的各级的出现概率P0(φ)、步骤S250A中计算的相互作用能φ的各级的出现概率P(φ)、步骤S300A中计算的∫Δν(φ)P(φ)dφ和数学式(1)来计算ΔG4a’，将所计算的ΔG4a’乘以负号，计算了ΔG4a。
[数学式12]
[式中，R表示气体常数，T表示通过分子动力学模拟设定的绝对温度(300开尔文)。]
需要说明的是，在数学式(1)的第1项中，通过将相互作用能φ进行单纯平均来算出。
所计算的ΔG4a为4.18(kcal/mol)。
(2) ΔG3a的计算
关于反应式(9)所表示的变化，以构成配体L0和位于该配体L0的周围的水分子的变化前的原子集合体作为“原子集合体A”、以构成配体L1和位于该配体L1周围的水分子的变化后的原子集合体作为“原子集合体AB”， 根据数学式(1)，计算了变化前的原子集合体A的自由能和片段B的自由能的和、与变化后的原子集合体AB的自由能之差ΔG3a’。
[化学式13]
需要说明的是，ΔG3a是与反应式(10)所表示的变化相关的自由能变化量，因此在根据数学式(1)计算的ΔG3a’上乘以负号而获得的值相当于ΔG3a。
[化学式14]
在本实施例中，利用图1所示的计算装置，通过实行图2所示的步骤S110～S190、S210～S280、S300和S400所对应的步骤S110B～S190B、S210B～S280B、S300B和S400B，计算了ΔG3a’。在步骤S120B中实行了图3所示的步骤S121～S123所对应的步骤S121B～S123B，在步骤S130B中实行了图5所示的步骤S131b～S134b所对应的步骤S131B～S134B，在步骤S220B中实行了图6所示的步骤S221～S223所对应的步骤S221B～S223B。
[变化前的原子集合体A所包含的配体L0和变化后的原子集合体AB所包含的配体L1的相关数据的输入]
在实行步骤S110B和S210B之前的准备步骤中，输入部11将变化前的原子集合体A所包含的配体L0和变化后的原子集合体AB所包含的配体L1的相关数据输入到控制部12中。所输入的数据通过控制部12存储在存储部13中。所输入的数据是从ΔG4a的计算中使用的复合体L0R和L1R的相关数据中抽取与配体L0和配体L1相关的数据而获得的数据，坐标数据分别存在于立体结构文件3PTB_L0、3PTB_L1中。
[模拟条件的相关数据的输入]
在实行步骤S110B和S210B之前的准备步骤中，输入部11将模拟条件的相关数据输入到控制部12中。所输入的数据通过控制部12存储在存储部13中。所输入的数据中，除与蛋白R有关的数据以外，其它数据与ΔG4a的计算相同。
[步骤S110B：第1原子集合体模型的制作]
在步骤S110B中，控制部12利用可视化软件VMD读取存储在存储部13中的PDB文件(3PTB_L0)，然后，利用该VMD中安装的“solvate”功能，以构成配体L0的原子的重心坐标为原点，使产生815个水分子，制作立方体型的基本单元，制作了将变化前的原子集合体A进行模型化而获得的第1原子集合体模型。需要说明的是，关于各原子的势能参数，在构成配体L0和配体L1的各原子中使用了CGenFF，在构成水分子的各原子中使用了TIP3P。另外，配体L0和配体L1的各原子的势能参数利用MATCH程序进行了确定。
[步骤S120B：通过计算机模拟获取原子集合体A的原子]
[步骤S130B：获取原子集合体AB的坐标]
[步骤S140B：相互作用能φ的计算]
[步骤S150B：相互作用能φ的各级频率分布的制作]
[步骤S160B：相互作用能φ的各级的出现概率P0(φ)的计算]
[步骤S170B：相互作用能ε的计算]
[步骤S180B：相互作用能φ的各级中的、相互作用能ε的各级频率分布的制作]
[步骤S190B：相互作用能φ的各级中的、相互作用能ε的各级的出现概率P0’(ε；φ)的计算]
控制部12进行与步骤S120A～S190A相同的操作，实行了步骤S120B～S190B。
需要说明的是，关于分子动力学模拟的条件，除了未通过阶段升温来实施平衡化而是在能量最小化后实施温度300K/模拟时间100皮可秒的平衡化模拟、以及将平衡化后模拟时间设为1000皮可秒、并以10飞秒的间隔输出原子集合体A的坐标以外，进行了与步骤S120A相同的操作。
[步骤S210B：第2原子集合体模型的制作]
[步骤S220B：通过计算机模拟获取原子集合体AB的坐标]
[步骤S230B：相互作用能φ的计算]
[步骤S240B：相互作用能φ的各级频率分布的制作]
[步骤S250B：相互作用能φ的各级的出现概率P(φ)的计算]
[步骤S260B：相互作用能ε的计算]
[步骤S270B：相互作用能ε的各级频率分布的制作]
[步骤S280B：相互作用能ε的各级的出现概率P’(ε)的计算]
控制部12进行与步骤S210A～S280A相同的操作，实行了步骤S210B～S280B。
需要说明的是，关于分子动力学模拟的条件，除了未通过阶段升温实施平衡化而是在能量最小化后实施温度300K/模拟时间100皮可秒的平衡化模拟、以及将平衡化后模拟时间设为1000皮可秒、并以10飞秒的间隔输出原子集合体AB的坐标以外，进行了与步骤S220A相同的操作。
[步骤S300B：相互作用能ε所引起的自由能变化量∫Δν(φ)P(φ)dφ的计算]
[步骤S400B：ΔG3a的计算]
控制部12进行与步骤S300A～S400A相同的操作，实行了步骤S300B～S400B。
所计算的ΔG3a为1.59(kcal/mol)。
(3) ΔG4b的计算
关于反应式(11)所表示的变化，以构成复合体L0R和位于该复合体L0R周围的水分子的变化前的原子集合体作为“原子集合体A”、以构成复合体L2R和位于该复合体L2R周围的水分子的变化后的原子集合体作为“原子集合体AB”，计算了变化前的原子集合体A的自由能和片段B的自由能的和、与变化后的原子集合体AB的自由能之差ΔG4b。
[化学式15]
在本实施例中，利用图1所示的计算装置，通过实行图2所示的步骤S110～S190、S210～S280、S300和S400所对应的步骤S110C～S190C、S210C～S280C、S300C和S400C，计算了ΔG4b。在步骤S120C中实行了图3所示的步骤S121～S123所对应的步骤S121C～S123C，在步骤S130C中实行了图5所示的步骤S131b～S134b所对应的步骤S131C～S134C，在步骤S220C中实行了图6所示的步骤S221～S223所对应的步骤S221C～S223C。
配体L2的化学结构如下。
[化学式16]
构成配体L2的各原子的部分电荷如下。
[化学式17]
片段B的结构如下。
[化学式18]
步骤S110C～S190C、S210C～S280C、S300C和S400C除了将配体L1变更为配体L2以及变更片段B以外，按照与ΔG4a的计算中的步骤S110A～S190A、S210A～S280A、S300A和S400A相同的方式来实行。所计算的ΔG4b为-17.84(kcal/mol)。
(4) ΔG3b的计算
关于反应式(12)所表示的变化，以构成配体L0和位于该配体L0周围的水分子的变化前的原子集合体为“原子集合体A”、以构成配体L2和位于该配体L2周围的水分子的变化后的原子集合体为“原子集合体AB”，计算了变化前的原子集合体A的自由能和片段B的自由能的和、与变化后的原子集合体AB的自由能之差ΔG3b。
[化学式19]
在本实施例中，利用图1所示的计算装置，通过实行图2所示的步骤S110～S190、S210～S280、S300和S400所对应的步骤S110D～S190D、S210D～S280D、S300D和S400D，计算了ΔG3b。在步骤S120D中实行了图3所示的步骤S121～S123所对应的步骤S121D～S123D，在步骤S130D中实行了图5所示的步骤S131b～S134b所对应的步骤S131D～S134D，在步骤S220D中实行了图6所示的步骤S221～S223所对应的步骤S221D～S223D。
步骤S110D～S190D、S210D～S280D、S300D和S400D除了将配体L1变更为配体L2并变更片段B以外，按照与ΔG3a的计算中的步骤S110B～S190B、S210B～S280B、S300B和S400B相同的方式来实行。所计算的ΔG3b为-14.91(kcal/mol)。
(5) ΔΔG的计算
由下述各式计算的ΔΔG为-0.34(kcal/mol)。
公知文献(Journal of Chemical Theory and Computation, 8, 3686-3695(2012)和Jouranl of the American Chemical Society, 99, 2331-2336(1977))中记载的、通过实验求得的ΔΔG为-0.10(kcal/mol)，通过本发明的计算方法计算的ΔΔG与通过实验求得的ΔΔG之差在1kcal/mol以内，一致性良好。
[实施例2]
使用配体L3来代替配体L2，进行与实施例1相同的操作，计算了ΔΔG。
配体L3的化学结构如下。
[化学式20]
构成配体L3的各原子的部分电荷如下。
[化学式21]
片段B的结构如下。
[化学式22]
这里，如下所述，为方便起见，将点电荷(部分电荷；0.000)作为片段B的原子来处理，从而按照与实施例1相同的方式计算了ΔΔG。
[化学式23]
ΔG4a和ΔG3a采用了与实施例1相同的值。
所计算的ΔG4b为-8.12(kcal/mol)。
所计算的ΔG3b为-5.45(kcal/mol)。
所计算的ΔΔG为-0.08(kcal/mol)。
公知文献(Journal of Chemical Theory and Computation, 8, 3686-3695(2012)和Jouranl of the American Chemical Society, 99, 2331-2336(1977))中记载的、通过实验求得的ΔΔG为0.27(kcal/mol)，通过本发明的计算方法算出的ΔΔG与通过实验求得的ΔΔG之差在1kcal/mol以内，一致性良好。
[实施例3]
使用配体L4来配体L2，进行与实施例1相同的操作，计算了ΔΔG。
配体L4的化学结构如下。
[化学式24]
构成配体L4的各原子的部分电荷如下。
[化学式25]
片段B的结构如下。
[化学式26]
ΔG4a和ΔG3a采用了与实施例1相同的值。
所计算的ΔG4b为-23.00(kcal/mol)。
所计算的ΔG3b为-19.88(kcal/mol)。
所计算的ΔΔG为-0.53(kcal/mol)。
公知文献(Journal of Chemical Theory and Computation, 8, 3686-3695(2012)和Jouranl of the American Chemical Society, 99, 2331-2336(1977))中记载的、通过实验求得的ΔΔG为-0.40(kcal/mol)，通过本发明的计算方法算出的ΔΔG与通过实验求得的ΔΔG之差在1kcal/mol以内，一致性良好。
[实施例4]
使用配体L5来代替配体L2，进行与实施例1相同的操作，计算了ΔΔG。
配体L5的化学结构如下。
[化学式27]
构成配体L5的各原子的部分电荷如下。
[化学式28]
片段B的结构如下。
[化学式29]
这里，如下所述，为方便起见，将点电荷(部分电荷；0.000)作为片段B的原子来处理，从而按照与实施例1相同的方式计算了ΔΔG。
[化学式30]
ΔG4a和ΔG3a采用了与实施例1相同的值。
所计算的ΔG4b为1.78(kcal/mol)。
所计算的ΔG3b为3.55(kcal/mol)。
所计算的ΔΔG为0.82(kcal/mol)。
公知文献(Jouranl of the American Chemical Society, 125, 10570-10579(2003))中记载的、通过实验求得的ΔΔG为0.93(kcal/mol)，通过本发明的计算方法算出的ΔΔG与通过实验求得的ΔΔG之差在1kcal/mol以内，一致性良好。
[实施例5]
使用配体L6来代替配体L2，进行与实施例1相同的操作，计算了ΔΔG。
配体L6的化学结构如下。
[化学式31]
构成配体L6的各原子的部分电荷如下。
[化学式32]
片段B的结构如下。
[化学式33]
这里，如下所述，为方便起见，将点电荷(部分电荷；0.000)作为片段B的原子来处理，从而按照与实施例1相同的方式计算了ΔΔG。
[化学式34]
ΔG4a和ΔG3a采用了与实施例1相同的值。
所计算的ΔG4b为-3.79(kcal/mol)。
所计算的ΔG3b为-0.95(kcal/mol)。
所计算的ΔΔG为-0.25(kcal/mol)。
公知文献(Jouranl of the American Chemical Society, 125, 10570-10579(2003))中记载的通过实验求出的ΔΔG为0.25(kcal/mol)，通过本发明的计算方法算出的ΔΔG与通过实验求得的ΔΔG之差在1kcal/mol以内，一致性良好。
符号说明
1：计算装置
11：输入部
12：控制部
12A：第1原子集合体模型制作单元
12B：第1坐标获取单元
12C：第2坐标获取单元
12D：第1相互作用能φ频率分布制作单元
12E：第1相互作用能φ出现概率计算单元
12F：第1相互作用能ε频率分布制作单元
12G：第1相互作用能ε出现概率计算单元
12H：第2原子集合体模型制作单元
12I：第3坐标获取单元
12J：第2相互作用能φ频率分布制作单元
12K：第2相互作用能φ出现概率计算单元
12L：第2相互作用能ε频率分布制作单元
12M：第2相互作用能ε出现概率计算单元
12N：∫Δν(φ)P(φ)dφ计算单元
12O：ΔG计算单元
13：存储部
13A：变化前的原子集合体A的相关数据
13B：第1模拟条件的相关数据
13C：第1模拟程序
13D：变化后的原子集合体AB的相关数据
13E：第2模拟条件的相关数据
13F：第2模拟程序
14：输出部

Claims (13)

1.一种计算装置，其具备控制部，关于反应式(1)所表示的变化，上述控制部计算变化前的原子集合体A的自由能和片段B的自由能的和、与变化后的原子集合体AB的自由能之差ΔG，
式中，A表示由结构a构成或包含结构a而构成的原子集合体，B表示由1个以上的原子构成的片段，AB表示由原子集合体A和与该原子集合体A的结构a结合的片段B构成的原子集合体，
其中，控制部具备下述单元：
第1原子集合体模型制作单元，制作将变化前的原子集合体A进行模型化而获得的第1原子集合体模型；
第1坐标获取单元，通过针对由第1原子集合体模型制作单元制作的第1原子集合体模型的计算机模拟，获取第1～第i的状态F1～Fi的各状态下的原子集合体A的坐标，其中i为2以上的整数；
第2坐标获取单元，根据通过第1坐标获取单元获取的原子集合体A的坐标，获取通过片段B与状态F1～Fi的各状态下的原子集合体A结合而产生的原子集合体AB的坐标；
第1相互作用能φ频率分布制作单元，根据通过第2坐标获取单元获取的原子集合体AB的坐标，计算结构a和与该结构a结合的片段B的相互作用能φ，制作表示相互作用能φ的各级频率的频率分布；
第1相互作用能φ出现概率计算单元，根据通过第1相互作用能φ频率分布制作单元制作的频率分布，计算相互作用能φ的各级的出现概率P0(φ)；
第1相互作用能ε频率分布制作单元，根据通过第2坐标获取单元获取的原子集合体AB的坐标，计算通过从原子集合体AB中去除由结构a和与该结构a结合的片段B构成的结构aB而产生的原子集合体的一部分或全部与片段B的相互作用能ε，制作通过第1相互作用能φ频率分布制作单元制作的频率分布的相互作用能φ的各级中的、表示相互作用能ε的各级频率的频率分布；
第1相互作用能ε出现概率计算单元，根据通过第1相互作用能ε频率分布制作单元制作的频率分布，计算在相互作用能φ的各级中的、相互作用能ε的各级的出现概率P0’(ε；φ)；
第2原子集合体模型制作单元，制作将变化后的原子集合体AB进行模型化而获得的第2原子集合体模型；
第3坐标获取单元，通过针对由第2原子集合体模型制作单元制作的第2原子集合体模型的计算机模拟，获取第1～第j的状态G1～Gj的各状态下的原子集合体AB的坐标，其中j为2以上的整数；
第2相互作用能φ频率分布制作单元，根据通过第3坐标获取单元获取的原子集合体AB的坐标，计算结构a和与该结构a结合的片段B的相互作用能φ，制作表示相互作用能φ的各级频率的频率分布；
第2相互作用能φ出现概率计算单元，根据通过第2相互作用能φ频率分布制作单元制作的频率分布，计算相互作用能φ的各级的出现概率P(φ)；
第2相互作用能ε频率分布制作单元，根据通过第3坐标获取单元获取的原子集合体AB的坐标，计算通过从原子集合体AB中去除由结构a和与该结构a结合的片段B构成的结构aB而产生的原子集合体的一部分或全部与片段B的相互作用能ε，制作表示相互作用能ε的各级频率的频率分布；
第2相互作用能ε出现概率计算单元，根据通过第2相互作用能ε频率分布制作单元制作的频率分布，计算相互作用能ε的各级的出现概率P’(ε)；
∫Δν(φ)P(φ)dφ计算单元，根据通过第2相互作用能φ出现概率计算单元计算的P(φ)、通过第1相互作用能ε出现概率计算单元计算的P0’(ε；φ)、以及通过第2相互作用能ε出现概率计算单元计算的P’(ε)，计算相互作用能ε所引起的自由能变化量∫Δν(φ)P(φ)dφ，式中Δν(φ)表示相互作用能φ的各级中的、相互作用能ε所引起的自由能变化量；以及
ΔG计算单元，根据通过第1相互作用能φ出现概率计算单元计算的P0(φ)、通过第2相互作用能φ出现概率计算单元计算的P(φ)、通过∫Δν(φ)P(φ)dφ计算单元计算的∫Δν(φ)P(φ)dφ和数学式(1)，来计算ΔG：
[数学式1]
式中，R表示气体常数，T表示反应式(1)所表示的变化发生时的绝对温度。
2.权利要求1所述的计算装置，其特征在于，第2坐标获取单元进行如下事项：
制作将原子集合体C进行模型化而获得的第3原子集合体模型，上述原子集合体C由结构a和与该结构a结合的片段B构成或者包含结构a和与该结构a结合的片段B而构成；
通过针对所制作的第3原子集合体模型的计算机模拟，获取第1～第k的状态H1～Hk的各状态下的原子集合体C的坐标，其中上述k为2以上的整数；
选择由选自构成结构a的原子的1个或2个以上的原子构成的选择原子组，使选自状态H1～Hk的1或2个以上的状态下的原子集合体C的选择原子组的坐标相对于状态F1～Fi的各状态下的原子集合体A的选择原子组的坐标旋转和/或并进，从而在原子集合体A的选择原子组与原子集合体C的选择原子组之间制作对应的原子间距离的平方和达到最小的原子集合体C的坐标，根据所制作的原子集合体C的坐标，使选自状态H1～Hk的1个或2个以上的状态下的原子集合体C与原子集合体A重叠，根据原子集合体A的坐标和与原子集合体A重叠的原子集合体C的片段B的1个或2个以上的坐标，获取通过片段B与原子集合体A结合而产生的1个或2个以上的原子集合体AB的坐标。
3.权利要求1或2所述的计算装置，其特征在于，∫Δν(φ)P(φ)dφ计算单元根据通过第2相互作用能φ出现概率计算单元计算的P(φ)、通过第1相互作用能ε出现概率计算单元计算的P0’(ε；φ)、以及通过第2相互作用能ε出现概率计算单元计算的P’(ε)，利用能量表示法计算相互作用能ε所引起的自由能变化量∫Δν(φ)P(φ)dφ。
4.权利要求1～3中任一项所述的计算装置，其特征在于，片段B由包含作为虚拟原子的点电荷的原子构成，第2坐标获取单元将片段B的点电荷添加在构成原子集合体A的结构a的原子的电荷参数中。
5.一种计算方法，关于反应式(1)所表示的变化，上述计算方法利用计算机计算变化前的原子集合体A的自由能和片段B的自由能的和、与变化后的原子集合体AB的自由能之差ΔG，
式中，A表示由结构a构成或包含结构a而构成的原子集合体，B表示由1个以上的原子构成的片段，AB表示由原子集合体A和与该原子集合体A的结构a结合的片段B构成的原子集合体，
其中，计算机的控制部实行如下步骤：
第1原子集合体模型制作步骤，制作将变化前的原子集合体A进行模型化而获得的第1原子集合体模型；
第1坐标获取步骤，通过针对第1原子集合体模型制作步骤中制作的第1原子集合体模型的计算机模拟，获取第1～第i的状态F1～Fi的各状态下的原子集合体A的坐标，其中i为2以上的整数；、
第2坐标获取步骤，根据通过第1坐标获取步骤获取的原子集合体A的坐标，获取通过片段B与状态F1～Fi的各状态下的原子集合体A结合而产生的原子集合体AB的坐标；
第1相互作用能φ频率分布制作步骤，根据通过第2坐标获取步骤获取的原子集合体AB的坐标，计算结构a和与该结构a结合的片段B的相互作用能φ，制作表示相互作用能φ的各级频率的频率分布；
第1相互作用能φ出现概率计算步骤，根据通过第1相互作用能φ频率分布制作步骤制作的频率分布，计算相互作用能φ的各级的出现概率P0(φ)；
第1相互作用能ε频率分布制作步骤，根据通过第2坐标获取步骤获取的原子集合体AB的坐标，计算通过从原子集合体AB中去除由结构a和与该结构a结合的片段B构成的结构aB而产生的原子集合体的一部分或全部与片段B的相互作用能ε，制作通过第1相互作用能φ频率分布制作步骤制作的频率分布的相互作用能φ的各级中的、表示相互作用能ε的各级频率的频率分布；
第1相互作用能ε出现概率计算步骤，根据通过第1相互作用能ε频率分布制作步骤制作的频率分布，计算相互作用能φ的各级中的、相互作用能ε的各级的出现概率P0’(ε；φ)；
第2原子集合体模型制作步骤，制作将变化后的原子集合体AB进行模型化而获得的第2原子集合体模型；
第3坐标获取步骤，通过针对第2原子集合体模型制作步骤中制作的第2原子集合体模型的计算机模拟，获取第1～第j的状态G1～Gj的各状态下的原子集合体AB的坐标，其中j为2以上的整数；
第2相互作用能φ频率分布制作步骤，根据通过第3坐标获取步骤获取的原子集合体AB的坐标，计算结构a和与该结构a结合的片段B的相互作用能φ，制作表示相互作用能φ的各级频率的频率分布；
第2相互作用能φ出现概率计算步骤，根据通过第2相互作用能φ频率分布制作步骤制作的频率分布，计算相互作用能φ的各级的出现概率P(φ)；
第2相互作用能ε频率分布制作步骤，根据通过第3坐标获取步骤获取的原子集合体AB的坐标，计算通过从原子集合体AB中去除由结构a和与该结构a结合的片段B构成的结构aB而产生的原子集合体的一部分或全部与片段B的相互作用能ε，制作表示相互作用能ε的各级频率的频率分布；
第2相互作用能ε出现概率计算步骤，根据通过第2相互作用能ε频率分布制作步骤制作的频率分布，计算相互作用能ε的各级的出现概率P’(ε)；
∫Δν(φ)P(φ)dφ计算步骤，根据通过第2相互作用能φ出现概率计算步骤计算的P(φ)、通过第1相互作用能ε出现概率计算步骤计算的P0’(ε；φ)、以及通过第2相互作用能ε出现概率计算步骤计算的P’(ε)，计算相互作用能ε所引起的自由能变化量∫Δν(φ)P(φ)dφ，式中Δν(φ)表示相互作用能φ的各级中的、相互作用能ε所引起的自由能变化量；以及
ΔG计算步骤，根据通过第1相互作用能φ出现概率计算步骤计算的P0(φ)、通过第2相互作用能φ出现概率计算步骤计算的P(φ)、通过∫Δν(φ)P(φ)dφ计算步骤计算的∫Δν(φ)P(φ)dφ和数学式(1)，来计算ΔG；
[数学式2]
式中，R表示气体常数，T表示反应式(1)所表示的变化发生时的绝对温度。
6.权利要求5所述的计算方法，其特征在于，在第2坐标获取步骤中，计算机的控制部实行如下事项：
制作将原子集合体C进行模型化而获得的第3原子集合体模型，上述原子集合体C由结构a和与该结构a结合的片段B构成或者包含结构a和与该结构a结合的片段B而构成；
通过针对所制作的第3原子集合体模型的计算机模拟，获取第1～第k的状态H1～Hk的各状态下的原子集合体C的坐标，其中k为2以上的整数；
选择由选自构成结构a的原子的1个或2个以上的原子构成的选择原子组，使选自状态H1～Hk的1种或2种以上的状态下的原子集合体C的选择原子组的坐标相对于状态F1～Fi的各状态下的原子集合体A的选择原子组的坐标旋转和/或并进，从而在原子集合体A的选择原子组与原子集合体C的选择原子组之间制作对应的原子间距离的平方和达到最小的原子集合体C的坐标，根据所制作的原子集合体C的坐标，使选自状态H1～Hk的1种或2种以上的状态下的原子集合体C与原子集合体A重叠，根据原子集合体A的坐标和与原子集合体A重叠的原子集合体C的片段B的1个或2个以上的坐标，获取通过片段B与原子集合体A结合而产生的1个或2个以上的原子集合体AB的坐标。
7.权利要求5或6所述的计算方法，其特征在于，在∫Δν(φ)P(φ)dφ计算步骤中，计算机的控制部根据通过第2相互作用能φ出现概率计算步骤计算的P(φ)、通过第1相互作用能ε出现概率计算步骤计算的P0’(ε；φ)、以及通过第2相互作用能ε出现概率计算步骤计算的P’(ε)，利用能量表示法计算相互作用能ε所引起的自由能变化量∫Δν(φ)P(φ)dφ。
8.权利要求5～7中任一项所述的计算方法，其特征在于，片段B由包含作为虚拟原子的点电荷的原子构成，在第2坐标获取步骤中，计算机的控制部将片段B的点电荷添加在构成原子集合体A的结构a的原子的电荷参数中。
9.一种程序，关于反应式(1)所表示的变化，上述程序计算变化前的原子集合体A的自由能和片段B的自由能的和、与变化后的原子集合体AB的自由能之差ΔG，
式中，A表示由结构a构成或者包含结构a而构成的原子集合体，B表示由1个以上的原子构成的片段，AB表示由原子集合体A和与该原子集合体A的结构a结合的片段B构成的原子集合体，
其中，上述程序使计算机的控制部起到如下单元的作用：
第1原子集合体模型制作单元，制作将变化前的原子集合体A进行模型化而获得的第1原子集合体模型；
第1坐标获取单元，通过针对由第1原子集合体模型制作单元制作的第1原子集合体模型的计算机模拟，获取第1～第i的状态F1～Fi的各状态下的原子集合体A的坐标，其中i为2以上的整数；
第2坐标获取单元，根据通过第1坐标获取单元获取的原子集合体A的坐标，获取通过片段B与状态F1～Fi的各状态下的原子集合体A结合而产生的原子集合体AB的坐标；
第1相互作用能φ频率分布制作单元，根据通过第2坐标获取单元获取的原子集合体AB的坐标，计算结构a和与该结构a结合的片段B的相互作用能φ，制作表示相互作用能φ的各级频率的频率分布；
第1相互作用能φ出现概率计算单元，根据通过第1相互作用能φ频率分布制作单元制作的频率分布，计算相互作用能φ的各级的出现概率P0(φ)；
第1相互作用能ε频率分布制作单元，根据通过第2坐标获取单元获取的原子集合体AB的坐标，计算通过从原子集合体AB中去除由结构a和与该结构a结合的片段B构成的结构aB而产生的原子集合体的一部分或全部与片段B的相互作用能ε，制作通过第1相互作用能φ频率分布制作单元制作的频率分布的相互作用能φ的各级中的、表示相互作用能ε的各级频率的频率分布；
第1相互作用能ε出现概率计算单元，根据通过第1相互作用能ε频率分布制作单元制作的频率分布，计算在相互作用能φ的各级中的、相互作用能ε的各级的出现概率P0’(ε；φ)；
第2原子集合体模型制作单元，制作将变化后的原子集合体AB进行模型化而获得的第2原子集合体模型；
第3坐标获取单元，通过针对由第2原子集合体模型制作单元制作的第2原子集合体模型的计算机模拟，获取第1～第j的状态G1～Gj的各状态下的原子集合体AB的坐标，其中j为2以上的整数；
第2相互作用能φ频率分布制作单元，根据通过第3坐标获取单元获取的原子集合体AB的坐标，计算结构a和与该结构a结合的片段B的相互作用能φ，制作表示相互作用能φ的各级频率的频率分布；
第2相互作用能φ出现概率计算单元，根据通过第2相互作用能φ频率分布制作单元制作的频率分布，计算相互作用能φ的各级的出现概率P(φ)；
第2相互作用能ε频率分布制作单元，根据通过第3坐标获取单元获取的原子集合体AB的坐标，计算通过从原子集合体AB中去除由结构a和与该结构a结合的片段B构成的结构aB而产生的原子集合体的一部分或全部与片段B的相互作用能ε，制作表示相互作用能ε的各级频率的频率分布；
第2相互作用能ε出现概率计算单元，根据通过第2相互作用能ε频率分布制作单元制作的频率分布，计算相互作用能ε的各级的出现概率P’(ε)；
∫Δν(φ)P(φ)dφ计算单元，根据通过第2相互作用能φ出现概率计算单元计算的P(φ)、通过第1相互作用能ε出现概率计算单元计算的P0’(ε；φ)、以及通过第2相互作用能ε出现概率计算单元计算的P’(ε)，计算相互作用能ε所引起的自由能变化量∫Δν(φ)P(φ)dφ，式中Δν(φ)表示相互作用能φ的各级中的、相互作用能ε所引起的自由能变化量；以及
ΔG计算单元，根据通过第1相互作用能φ出现概率计算单元计算的P0(φ)、通过第2相互作用能φ出现概率计算单元计算的P(φ)、通过∫Δν(φ)P(φ)dφ计算单元计算的∫Δν(φ)P(φ)dφ和数学式(1)，来计算ΔG：
[数学式3]
式中，R表示气体常数，T表示反应式(1)所表示的变化发生时的绝对温度。
10.权利要求9所述的程序，其特征在于，第2坐标获取单元实行如下事项：
制作将原子集合体C进行模型化而获得的第3原子集合体模型，上述原子集合体C由结构a和与该结构a结合的片段B构成或者包含结构a和与该结构a结合的片段B而构成；
通过针对所制作的第3原子集合体模型的计算机模拟，获取第1～第k的状态H1～Hk的各状态下的原子集合体C的坐标，其中k为2以上的整数；
选择由选自构成结构a的原子的1个或2个以上的原子构成的选择原子组，使选自状态H1～Hk的1种或2种以上的状态下的原子集合体C的选择原子组的坐标相对于状态F1～Fi的各状态下的原子集合体A的选择原子组的坐标旋转和/或并进，从而在原子集合体A的选择原子组与原子集合体C的选择原子组之间制作对应的原子间距离的平方和达到最小的原子集合体C的坐标，根据所制作的原子集合体C的坐标，使选自状态H1～Hk的1种或2种以上的状态下的原子集合体C与原子集合体A重叠，根据原子集合体A的坐标和与原子集合体A重叠的原子集合体C的片段B的1个或2个以上的坐标，获取通过片段B与原子集合体A结合而产生的1个或2个以上的原子集合体AB的坐标。
11.权利要求9或10所述的程序，其特征在于，Δν(φ)P(φ)dφ计算单元根据通过第2相互作用能φ出现概率计算单元计算的P(φ)、通过第1相互作用能ε出现概率计算单元计算的P0’(ε；φ)、以及通过第2相互作用能ε出现概率计算单元计算的P’(ε)，利用能量表示法计算相互作用能ε所引起的自由能变化量∫Δν(φ)P(φ)dφ。
12.权利要求9～11中任一项所述的程序，其特征在于，片段B由包含作为虚拟原子的点电荷的原子构成，第2坐标获取单元将片段B的点电荷添加在构成原子集合体A的结构a的原子的电荷参数中。
13.一种计算机可读取的记录介质，其记录有权利要求9～12中任一项所述的程序。
CN201580052788.3A 2014-09-30 2015-09-30 自由能计算装置、方法、程序、以及记录有该程序的记录介质 Active CN107004064B (zh)

Applications Claiming Priority (3)

Application Number Priority Date Filing Date Title
JP2014202658 2014-09-30
JP2014-202658 2014-09-30
PCT/JP2015/077847 WO2016052662A1 (ja) 2014-09-30 2015-09-30 自由エネルギー計算装置、方法、プログラム、並びに該プログラムを記録した記録媒体

Publications (2)

Publication Number Publication Date
CN107004064A true CN107004064A (zh) 2017-08-01
CN107004064B CN107004064B (zh) 2020-09-29

Family

ID=55630685

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201580052788.3A Active CN107004064B (zh) 2014-09-30 2015-09-30 自由能计算装置、方法、程序、以及记录有该程序的记录介质

Country Status (6)

Country Link
US (1) US11126761B2 (zh)
EP (1) EP3203401A4 (zh)
JP (1) JP6565102B2 (zh)
CN (1) CN107004064B (zh)
CA (1) CA2962730A1 (zh)
WO (1) WO2016052662A1 (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2023102870A1 (zh) * 2021-12-10 2023-06-15 上海智药科技有限公司 分子自由能计算、稳定性分析方法、装置、设备及存储介质

Families Citing this family (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20150178442A1 (en) 2013-12-23 2015-06-25 Schrodinger, Inc. Methods and systems for calculating free energy differences using a modified bond stretch potential
US10726946B2 (en) * 2017-08-22 2020-07-28 Schrödinger, Inc. Methods and systems for calculating free energy differences using an alchemical restraint potential
CN114187971B (zh) * 2021-12-10 2023-03-28 上海智药科技有限公司 分子自由能计算、稳定性分析方法、装置、设备及存储介质
JP2023110211A (ja) * 2022-01-28 2023-08-09 国立大学法人北海道大学 反応経路探索プログラム、反応経路探索システム及び反応経路探索方法

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102272762A (zh) * 2010-03-24 2011-12-07 松下电器产业株式会社 相互作用力变化预测装置和相互作用力变化预测方法
CN103049677A (zh) * 2011-10-17 2013-04-17 同济大学 一种处理小分子和蛋白质相互作用的计算方法
WO2013142630A1 (en) * 2012-03-20 2013-09-26 University Of Maryland, Baltimore Site-specific fragment identification guided by single-step free energy perturbation calculations
EP2782033A1 (en) * 2013-03-22 2014-09-24 Fujitsu Limited Calculation method of binding free energy, calculation device of binding free energy, program, screening method of compound

Family Cites Families (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5854992A (en) * 1996-09-26 1998-12-29 President And Fellows Of Harvard College System and method for structure-based drug design that includes accurate prediction of binding free energy
JP5673245B2 (ja) 2011-03-14 2015-02-18 富士通株式会社 自由エネルギー差予測方法及びシミュレーション装置
WO2012174208A1 (en) 2011-06-14 2012-12-20 Florida State University Research Foundation, Inc. Methods and apparatus for double-integration orthogonal space tempering

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102272762A (zh) * 2010-03-24 2011-12-07 松下电器产业株式会社 相互作用力变化预测装置和相互作用力变化预测方法
CN103049677A (zh) * 2011-10-17 2013-04-17 同济大学 一种处理小分子和蛋白质相互作用的计算方法
WO2013142630A1 (en) * 2012-03-20 2013-09-26 University Of Maryland, Baltimore Site-specific fragment identification guided by single-step free energy perturbation calculations
EP2782033A1 (en) * 2013-03-22 2014-09-24 Fujitsu Limited Calculation method of binding free energy, calculation device of binding free energy, program, screening method of compound

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2023102870A1 (zh) * 2021-12-10 2023-06-15 上海智药科技有限公司 分子自由能计算、稳定性分析方法、装置、设备及存储介质

Also Published As

Publication number Publication date
EP3203401A4 (en) 2018-08-08
EP3203401A1 (en) 2017-08-09
CA2962730A1 (en) 2016-04-07
JPWO2016052662A1 (ja) 2017-07-27
CN107004064B (zh) 2020-09-29
WO2016052662A1 (ja) 2016-04-07
JP6565102B2 (ja) 2019-08-28
US11126761B2 (en) 2021-09-21
US20170220717A1 (en) 2017-08-03

Similar Documents

Publication Publication Date Title
CN107004064A (zh) 自由能计算装置、方法、程序、以及记录有该程序的记录介质
Patkowski Recent developments in symmetry‐adapted perturbation theory
Tong et al. Accelerating CALYPSO structure prediction by data-driven learning of a potential energy surface
Del Ben et al. Electron correlation in the condensed phase from a resolution of identity approach based on the Gaussian and plane waves scheme
Marti et al. Complete-graph tensor network states: a new fermionic wave function ansatz for molecules
Deckard et al. Preliminary studies on the in silico evolution of biochemical networks
Nguyen-Truong et al. An implementation of the Levenberg–Marquardt algorithm for simultaneous-energy-gradient fitting using two-layer feed-forward neural networks
Bespalova et al. Hamiltonian operator approximation for energy measurement and ground-state preparation
CN109885917B (zh) 一种并行分子动力学模拟方法及系统
Goez et al. Embedding methods in quantum chemistry
Reith et al. A modern workflow for force-field development–Bridging quantum mechanics and atomistic computational models
Jaśkowiec et al. High‐order symmetric cubature rules for tetrahedra and pyramids
Backhouse et al. Scalable and predictive spectra of correlated molecules with moment truncated iterated perturbation theory
CN106503483A (zh) 基于模块化因子图的骨髓瘤信号通路机制确认方法
Chen et al. Mapping environment-specific quantitative trait loci
Putz Quantum theory: density, condensation, and bonding
Kobayashi et al. Divide-and-conquer approaches to quantum chemistry: Theory and implementation
Kossoski et al. Seniority and hierarchy configuration interaction for radicals and excited states
Heidari et al. Concurrent coupling of realistic and ideal models of liquids and solids in Hamiltonian adaptive resolution simulations
CN105808508B (zh) 一种求解不确定热传导问题的随机正交展开方法
Morita et al. IMiCMO: a new integrated ab initio multicenter molecular orbitals method for molecular dynamics calculations in solvent cluster systems
Capelo Using species abundance and phylogeny conjointly to approach vegetation classification: A case study on Macaronesia’s woody vegetation
Kanso et al. Reducing a model of sugar metabolism in peach to catch different patterns among genotypes
Sternberg et al. NOON—a non-orthogonal localised orbital order-N method
Wildová et al. The contrasting roles of growth traits and architectural traits in diversity maintenance in clonal plant communities

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant