空间站广场

论文

Notebooks

比赛

课程

Apps

我的主页

我的Notebooks

我的论文库

我的足迹

我的工作空间

任务

节点

文件

数据集

镜像

项目

数据库

公开

使用 Poltype2 获得Amoeba力场参数

分子动力学

力场

拟合

分子动力学力场拟合

hanyanbo

发布于 2023-08-22

推荐镜像 :Basic Image:bohrium-notebook:2023-04-07

推荐机型 :c2_m4_cpu

数据集

输入文件(v1)

使用 Poltype2 获得Amoeba力场参数

©️ Copyright 2023 @ Authors
作者： 韩言博 📨
日期：2023-08-21
共享协议：本作品采用知识共享署名-非商业性使用-相同方式共享 4.0 国际许可协议进行许可。
快速开始：点击上方的 开始连接 按钮，选择 bohrium-notebook 系列镜像和任意配置机型即可开始。（notebook运行镜像 Python=3.8 为佳）

代码

文本

推荐前置知识

分子动力学 | 力场

dflow

代码

文本

简要介绍：AMOEBA力场

AMOEBA力场的形式

AMOEBA力场（全称Atomic Multipole Optimized Energetics for Biomolecular Applications），是由 Washington University (St. Louis) 的 Jay Ponder领衔开发的包含了对分子体系多极能量矫正的分子力场。与其他力场相比，AMOEBA力场最大的不同是向力场中引入了对电相互作用的更细致描述，如下式所述

$E = = E_{intra} E_{bond} + E_{angle} + E_{bond - angle - cross} + E_{out - of - plane} + E_{torsion} + + E_{inter} E_{vdW} + E_{ele}^{permanent} + E_{ele}^{induced}$

该部分的能量项用于更好地计算分子间因为极化作用产生的复杂静电相互作用，对于溶剂体系中分子之间存在方向性的极化作用要有更细致的描述。

代码

文本

提示：本作品会调用 Tinker 等采用TINKER LICENSE的工具，您应当自行阅读该协议，并仅在该协议的限制要求下使用。

代码

文本

参数拟合的说明

AMOEBA力场中的键参（bond、angle等项）和经典力场的参数基本可以对应，只是其形式上引入了更多高阶项，这里不赘述；描述其分子间相互作用的vdw也基本类似，而描述静电相互作用的多极部分可以通过 TinkerTools/poltype2 工具来实现拟合。其拟合标的是分子的表面静电势能面（Electrostatic Potential Surfaces），也是分子对外界电荷的响应，这一数值可以使用量子化学计算软件（例如Psi4）计算得到。通过拟合ESP，AMOEBA力场中的原子多极将具有对分子间极化作用的细致描述。

Poltype2中是配合Tinker的poledit来实现这个过程的，其拟合的详细步骤是：

读取sdf或mol格式的文件，并做检查
筛选出需要进行多极拟合的原子，并匹配原子类型，产生必要的构象
- 筛选过大的分子和环状分子，减少其构象的生成
- 计算对称性
- 在对称性基础上提取出原子类型，并生成多极结构的局域坐标
- 使用rdkit产生构象
- 考虑必要的电离结构和异构体（多见于变pH环境）
检查原子是否可以使用Psi4进行计算，除非有I、Br等原子，均使用MP2基组计算
- Psi4计算几何优化任务
- 并检查计算结果的拓扑没有改变
对照力场库，搜索键参、vdw等力场参数
通过GDMA计算生成多极分解初猜
高精度的量子化学计算生成ESP数据
调用poledit拟合多极矩，使之尽可能符合静电势
输出Tinker的key参数文件

代码

文本

正如其名，AMOEBA力场相对成熟的更多地是在生物分子的应用场景，因此poltype2核心支持的范围除典型的碳氢和氮氢键、氧氢和个别自由基外，是以下键类型和形式电荷：

引自Poltype2/PoltypeModules/poltype.py CheckInputCharge()函数（https://github.com/TinkerTools/poltype2/blob/master/PoltypeModules/poltype.py）

            atomicnumtoformalchg={1:{2:1},5:{4:1},6:{3:-1},7:{2:-1,4:1},8:{1:-1,3:1},15:{4:1},16:{1:-1,3:1,5:-1},17:{0:-1,4:3},9:{0:-1},35:{0:-1},53:{0:-1}}

即氢键、四配位硼、sp2碳、胺和铵、羟基、𨦡/水合质子/质子化羟基、膦、巯基和硫酰、酰氯和卤素氟氯溴碘的离子(不成键)。

代码

文本

拟合小分子的Amoeba力场

poltype2封装了大部分的工具，本文将参数和环境配置的部分进一步封装方便使用，只需要准备结构文件即可，这里以DMC为例说明。

Tips: 如果你希望调整输入参数，可参阅poltype2文档，修改poltype.ini文件的写入（https://github.com/TinkerTools/poltype2）对酯类醚类小分子，大部分情况下本文使用的参数足够适用。

代码

文本

拟合工作流

我们使用dflow构建来运行拟合，使用FittingOP来执行poltype的工作命令。该OP预期的输出为key_file、拟合日志文件和部分参数的检查参考图。

代码

文本

[1]

!pip install -U pydflow lbg

Looking in indexes: https://pypi.tuna.tsinghua.edu.cn/simple
Collecting pydflow
  Downloading https://pypi.tuna.tsinghua.edu.cn/packages/49/af/a466ad5cb54cee614c59a37b7cfca43b2be7f412c0d16af3647a49264ba4/pydflow-1.7.67-py3-none-any.whl (132 kB)
     ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 132.7/132.7 kB 5.3 MB/s eta 0:00:00
Requirement already satisfied: ipython in /opt/conda/lib/python3.8/site-packages (8.10.0)
Collecting ipython
  Downloading https://pypi.tuna.tsinghua.edu.cn/packages/18/39/d5911ef659dbae6398dfcb7ac9e266a8b886ed04f338ecf04b35f934ef87/ipython-8.12.2-py3-none-any.whl (797 kB)
     ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 797.8/797.8 kB 16.7 MB/s eta 0:00:00a 0:00:01
Collecting jsonpickle
  Downloading https://pypi.tuna.tsinghua.edu.cn/packages/d3/25/6e0a450430b7aa194b0f515f64820fc619314faa289458b7dfca4a026114/jsonpickle-3.0.2-py3-none-any.whl (40 kB)
     ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 40.7/40.7 kB 13.8 MB/s eta 0:00:00
Requirement already satisfied: certifi in /opt/conda/lib/python3.8/site-packages (from pydflow) (2022.12.7)
Requirement already satisfied: psutil in /opt/conda/lib/python3.8/site-packages (from pydflow) (5.9.0)
Requirement already satisfied: urllib3 in /opt/conda/lib/python3.8/site-packages (from pydflow) (1.26.14)
Requirement already satisfied: six in /opt/conda/lib/python3.8/site-packages (from pydflow) (1.16.0)
Collecting typeguard==2.13.3
  Downloading https://pypi.tuna.tsinghua.edu.cn/packages/9a/bb/d43e5c75054e53efce310e79d63df0ac3f25e34c926be5dffb7d283fb2a8/typeguard-2.13.3-py3-none-any.whl (17 kB)
Requirement already satisfied: pyyaml in /opt/conda/lib/python3.8/site-packages (from pydflow) (6.0)
Requirement already satisfied: requests in /opt/conda/lib/python3.8/site-packages (from pydflow) (2.28.2)
Collecting cloudpickle==2.2.0
  Downloading https://pypi.tuna.tsinghua.edu.cn/packages/cf/26/cd6c4177273ee35f7a31245893489c68bc340988f12ca315b392f1f18a93/cloudpickle-2.2.0-py3-none-any.whl (25 kB)
Requirement already satisfied: tqdm in /opt/conda/lib/python3.8/site-packages (from pydflow) (4.64.1)
Collecting minio
  Downloading https://pypi.tuna.tsinghua.edu.cn/packages/9a/ae/c20305102e32f078ec5e77a5010378318307a2de213f61be6c5facb56f0f/minio-7.1.16-py3-none-any.whl (77 kB)
     ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 77.9/77.9 kB 25.7 MB/s eta 0:00:00
Requirement already satisfied: python-dateutil in /opt/conda/lib/python3.8/site-packages (from pydflow) (2.8.2)
Collecting argo-workflows==5.0.0
  Downloading https://pypi.tuna.tsinghua.edu.cn/packages/b1/6a/8f13d5124b111e8e054594d23782ea9c5dadda0517d1dd9ad08c7c055732/argo_workflows-5.0.0-py3-none-any.whl (452 kB)
     ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 452.5/452.5 kB 29.0 MB/s eta 0:00:00
Collecting kubernetes
  Downloading https://pypi.tuna.tsinghua.edu.cn/packages/99/89/3ab0cb3069f49ae2eaf73f884c82164f18f70fcc598e0312edea71614ad7/kubernetes-27.2.0-py2.py3-none-any.whl (1.5 MB)
     ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 1.5/1.5 MB 33.1 MB/s eta 0:00:00a 0:00:01
Requirement already satisfied: stack-data in /opt/conda/lib/python3.8/site-packages (from ipython) (0.2.0)
Requirement already satisfied: pickleshare in /opt/conda/lib/python3.8/site-packages (from ipython) (0.7.5)
Requirement already satisfied: prompt-toolkit!=3.0.37,<3.1.0,>=3.0.30 in /opt/conda/lib/python3.8/site-packages (from ipython) (3.0.36)
Requirement already satisfied: matplotlib-inline in /opt/conda/lib/python3.8/site-packages (from ipython) (0.1.6)
Requirement already satisfied: pygments>=2.4.0 in /opt/conda/lib/python3.8/site-packages (from ipython) (2.11.2)
Requirement already satisfied: jedi>=0.16 in /opt/conda/lib/python3.8/site-packages (from ipython) (0.18.1)
Requirement already satisfied: typing-extensions in /opt/conda/lib/python3.8/site-packages (from ipython) (4.5.0)
Requirement already satisfied: decorator in /opt/conda/lib/python3.8/site-packages (from ipython) (5.1.1)
Requirement already satisfied: pexpect>4.3 in /opt/conda/lib/python3.8/site-packages (from ipython) (4.8.0)
Requirement already satisfied: traitlets>=5 in /opt/conda/lib/python3.8/site-packages (from ipython) (5.7.1)
Requirement already satisfied: backcall in /opt/conda/lib/python3.8/site-packages (from ipython) (0.2.0)
Requirement already satisfied: parso<0.9.0,>=0.8.0 in /opt/conda/lib/python3.8/site-packages (from jedi>=0.16->ipython) (0.8.3)
Requirement already satisfied: ptyprocess>=0.5 in /opt/conda/lib/python3.8/site-packages (from pexpect>4.3->ipython) (0.7.0)
Requirement already satisfied: wcwidth in /opt/conda/lib/python3.8/site-packages (from prompt-toolkit!=3.0.37,<3.1.0,>=3.0.30->ipython) (0.2.5)
Requirement already satisfied: google-auth>=1.0.1 in /opt/conda/lib/python3.8/site-packages (from kubernetes->pydflow) (2.16.1)
Requirement already satisfied: websocket-client!=0.40.0,!=0.41.*,!=0.42.*,>=0.32.0 in /opt/conda/lib/python3.8/site-packages (from kubernetes->pydflow) (1.5.1)
Requirement already satisfied: oauthlib>=3.2.2 in /opt/conda/lib/python3.8/site-packages (from kubernetes->pydflow) (3.2.2)
Requirement already satisfied: requests-oauthlib in /opt/conda/lib/python3.8/site-packages (from kubernetes->pydflow) (1.3.1)
Requirement already satisfied: idna<4,>=2.5 in /opt/conda/lib/python3.8/site-packages (from requests->pydflow) (3.4)
Requirement already satisfied: charset-normalizer<4,>=2 in /opt/conda/lib/python3.8/site-packages (from requests->pydflow) (3.0.1)
Requirement already satisfied: asttokens in /opt/conda/lib/python3.8/site-packages (from stack-data->ipython) (2.0.5)
Requirement already satisfied: pure-eval in /opt/conda/lib/python3.8/site-packages (from stack-data->ipython) (0.2.2)
Requirement already satisfied: executing in /opt/conda/lib/python3.8/site-packages (from stack-data->ipython) (0.8.3)
Requirement already satisfied: cachetools<6.0,>=2.0.0 in /opt/conda/lib/python3.8/site-packages (from google-auth>=1.0.1->kubernetes->pydflow) (5.3.0)
Requirement already satisfied: pyasn1-modules>=0.2.1 in /opt/conda/lib/python3.8/site-packages (from google-auth>=1.0.1->kubernetes->pydflow) (0.2.8)
Requirement already satisfied: rsa<5,>=3.1.4 in /opt/conda/lib/python3.8/site-packages (from google-auth>=1.0.1->kubernetes->pydflow) (4.9)
Requirement already satisfied: pyasn1<0.5.0,>=0.4.6 in /opt/conda/lib/python3.8/site-packages (from pyasn1-modules>=0.2.1->google-auth>=1.0.1->kubernetes->pydflow) (0.4.8)
Installing collected packages: argo-workflows, typeguard, minio, jsonpickle, cloudpickle, kubernetes, ipython, pydflow
  Attempting uninstall: typeguard
    Found existing installation: typeguard 2.7.1
    Uninstalling typeguard-2.7.1:
      Successfully uninstalled typeguard-2.7.1
  Attempting uninstall: cloudpickle
    Found existing installation: cloudpickle 2.2.1
    Uninstalling cloudpickle-2.2.1:
      Successfully uninstalled cloudpickle-2.2.1
  Attempting uninstall: ipython
    Found existing installation: ipython 8.10.0
    Uninstalling ipython-8.10.0:
      Successfully uninstalled ipython-8.10.0
ERROR: pip's dependency resolver does not currently take into account all the packages that are installed. This behaviour is the source of the following dependency conflicts.
ipywidgets 7.6.5 requires widgetsnbextension~=3.5.0, but you have widgetsnbextension 3.6.2 which is incompatible.
Successfully installed argo-workflows-5.0.0 cloudpickle-2.2.0 ipython-8.12.2 jsonpickle-3.0.2 kubernetes-27.2.0 minio-7.1.16 pydflow-1.7.67 typeguard-2.13.3
WARNING: Running pip as the 'root' user can result in broken permissions and conflicting behaviour with the system package manager. It is recommended to use a virtual environment instead: https://pip.pypa.io/warnings/venv

代码

文本

[6]

import shutil

from pathlib import Path

import getpass

from dflow import upload_artifact, Step, Workflow

from dflow.python import OP, OPIO, Artifact, OPIOSign, Parameter

from dflow.python import PythonOPTemplate

POLTYPE_COMMAND='''source /root/.bashrc

export myusername=`whoami`

conda activate amoebamdpoltype

export GDMADIR=/opt

export PATH=/opt:$PATH

export PSI_SCRATCH=/scratch/$myusername/

export PATH=/root/tinker/bin:$PATH

mkdir /scratch

mkdir /scratch/root

python /root/poltype2/PoltypeModules/poltype.py'''

class FittingOP(OP):

@classmethod

def get_input_sign(cls):

return OPIOSign(

{

'sdf_or_mol2_file': Artifact(Path),

'run_script': Parameter(str, default=POLTYPE_COMMAND)

}

)

@classmethod

def get_output_sign(cls):

return OPIOSign({

'key_file': Artifact(Path),

'log_file': Artifact(Path),

'plot_file': Artifact(Path),

})

@OP.exec_sign_check

def execute(self, op_in: OPIO) -> OPIO:

op_in['sdf_or_mol2_file']: Path

# copy sdf_or_mol2_file to current dir

shutil.copy(Path(op_in['sdf_or_mol2_file']), Path(op_in['sdf_or_mol2_file']).name)

# create ini file

Path('poltype.ini').write_text('''structure={}

addhydrogentocharged=False\n

generateextendedconf=True

'''.format(Path(op_in['sdf_or_mol2_file']).name))

# create bash script

run_script = op_in['run_script']

Path('run.sh').write_text(str(run_script))

# priveleges

import os

os.system('chmod +x run.sh')

# run script

os.system('bash run.sh >tmp_log')

keyfilepath = Path('final.key').rename(

Path('final.key').with_name(

Path(op_in['sdf_or_mol2_file']).with_suffix('.key').name

)

return OPIO({

'key_file': keyfilepath,

'log_file': Path('tmp_log'),

'plot_file': Path('OPENME'),

})

代码

文本

配置好dflow的运行环境，需要使用你的Bohrium账号和密码：

代码

文本

[2]

from dflow import config, s3_config

from dflow.plugins import bohrium

from dflow.plugins.bohrium import TiefblueClient

from dflow.plugins.dispatcher import DispatcherExecutor

username = getpass.getpass('Bohrium账号：')

password = getpass.getpass('Bohrium密码:')

project_id = getpass.getpass('Bohrium项目ID：')

config["host"] = "https://workflows.deepmodeling.com"

config["k8s_api_server"] = "https://workflows.deepmodeling.com"

bohrium.config["username"] = username

bohrium.config["password"] = password

bohrium.config["project_id"] = project_id

s3_config["repo_key"] = "oss-bohrium"

s3_config["storage_client"] = TiefblueClient()

fitting_bohrium_config = {

"dispatcher": {

"batch_type": "Bohrium",

"context_type": "Bohrium",

"input_data": {

"job_type": "container",

"platform": "ali",

"scass_type": "c16_m32_cpu",

}

Bohrium账号： ········
Bohrium密码: ········
Bohrium项目ID： ········

代码

文本

我们需要一个描述分子结构和拓扑信息的sdf文件作为输入，sdf文件可以通过以下代码写入，也可直接使用你自己的sdf文件。

sdf文件的格式是：

可空的分子名称
可空的注释
<空行>
原子数 键的个数
坐标信息起始行
...
坐标信息最后一行
键信息起始行
...
键键信息最后一行
M END

代码

文本

[4]

import logging

from pathlib import Path

sdf_content = """

12 11 0 0 0 0 0 0 0 0999 V2000

-0.3288 0.6210 -0.2919 C 0 0 0 0 0 0 0 0 0 0 0 0

0.1644 1.7089 -0.3122 O 0 0 0 0 0 0 0 0 0 0 0 0

0.2804 -0.4266 0.2499 O 0 0 0 0 0 0 0 0 0 0 0 0

-1.5387 0.3137 -0.7927 O 0 0 0 0 0 0 0 0 0 0 0 0

-2.2612 1.3720 -1.4078 C 0 0 0 0 0 0 0 0 0 0 0 0

-3.3149 1.1959 -1.2017 H 0 0 0 0 0 0 0 0 0 0 0 0

-1.9407 2.3313 -0.9988 H 0 0 0 0 0 0 0 0 0 0 0 0

-2.0818 1.3584 -2.4843 H 0 0 0 0 0 0 0 0 0 0 0 0

-0.3671 -1.6910 0.2448 C 0 0 0 0 0 0 0 0 0 0 0 0

-0.5639 -2.0241 -0.7757 H 0 0 0 0 0 0 0 0 0 0 0 0

0.3254 -2.3719 0.7346 H 0 0 0 0 0 0 0 0 0 0 0 0

-1.3076 -1.6493 0.7967 H 0 0 0 0 0 0 0 0 0 0 0 0

1 3 1 0 0 0 0

2 1 2 0 0 0 0

4 1 1 0 0 0 0

5 6 1 0 0 0 0

5 7 1 0 0 0 0

5 4 1 0 0 0 0

8 5 1 0 0 0 0

9 3 1 0 0 0 0

9 11 1 0 0 0 0

9 12 1 0 0 0 0

10 9 1 0 0 0 0

M END

$$$$

"""

sdf_or_mol2_file_path = Path("DMC.sdf") # 可以修改为真实的sdf文件地址

# for directly user to avoid file overwrite

if sdf_or_mol2_file_path.exists():

logging.warning(f'{sdf_or_mol2_file_path} exists!')

overwrite = '?'

while overwrite not in ('y', '', 'N', 'n'):

overwrite = input('override?(y/N):')

if overwrite == 'y':

sdf_or_mol2_file_path.write_text(sdf_content)

logging.warning('File override.')

elif overwrite in ('N','','n'):

logging.warning('Pass override.')

continue

else:

logging.warning('Wrong choice {}!'.format(overwrite))

else:

sdf_or_mol2_file_path.write_text(sdf_content)

logging.warning('File written.')

WARNING:root:Path("DMC.sdf") exists!
override?(y/N): 
WARNING:root:Pass override.

代码

文本

[ ]

sdf_or_mol2_file = upload_artifact(sdf_or_mol2_file_path)

fitting_step = Step(

name='fitting-step',

template=PythonOPTemplate(

FittingOP,

command=[

# 实际上你也可以在这里 hack 形式电荷，但多数情况下用处不大

# "sed -i \\\"s/atomicnumtoformalchg\\\\={1:{2:1},5:{4:1},6:{3:-1},7:{2:-1,4:1},8:{1:-1,3:1},15:{4:1},16:{1:-1,3:1,5:-1},17:{0:-1,4:3},9:{0:-1},35:{0:-1},53:{0:-1}}/atomicnumtoformalchg\\\\={1:{2:1},5:{4:1},6:{3:-1},7:{2:-1,4:1},8:{1:-1,3:1},15:{4:1,6:-1},16:{1:-1,3:1,5:-1},17:{0:-1,4:3},9:{0:-1},35:{0:-1},53:{0:-1}}/\\\" /root/poltype2/PoltypeModules/poltype.py && sed -i \\\"s/Chem.SanitizeMol(rdkitmol)//g\\\" /root/poltype2/PoltypeModules/symmetry.py && conda run -n amoebamdpoltype python"

"conda run -n amoebamdpoltype python"

image='registry.dp.tech/dptech/prod-13211/poltype2-cpu:159f352c-cpu-xtb',

executor=DispatcherExecutor(

machine_dict={

"batch_type": fitting_bohrium_config["dispatcher"].get("batch_type"),

"context_type": fitting_bohrium_config["dispatcher"].get("context_type"),

"remote_profile": {"input_data": fitting_bohrium_config["dispatcher"].get("input_data")}

image_pull_policy = "IfNotPresent"),

artifacts={

'sdf_or_mol2_file': sdf_or_mol2_file,

)

wf = Workflow(name='poltype2-{}'.format(str(sdf_or_mol2_file_path.stem[:5]).lower()))

wf.add(fitting_step)

wf.submit();

logging.info("Please wait the workflow to finish...")

wf.wait()

logging.info("Workflow done.")

代码

文本

如果拟合正常结束，那么你将可以下载到输出文件，分别是tinker可用的力场参数文件.key、任务日志文件和对力场拟合效果的benchmark图片。

代码

文本

[ ]

from dflow import download_artifact

download_artifact(wf.query_step(name="fitting-step")[0].outputs.artifacts["key_file"])

download_artifact(wf.query_step(name="fitting-step")[0].outputs.artifacts["log_file"])

download_artifact(wf.query_step(name="fitting-step")[0].outputs.artifacts["plot_file"])

代码

文本

下载好文件后，可以看到key文件中的各项力场参数，以及他们的来源

也能看到，除了最后的multipole参数，都是从现有的力场库中匹配得到的

代码

文本

[18]

print_key = input(".key file will be long, print it?(Y/N):")

if print_key in ('Y','y'):

print(Path(sdf_or_mol2_file_path.with_suffix('.key')).read_text())

.key file will be long, print it?(Y/N): Y
parameters /root/poltype2/ParameterFiles/amoebabio18_header.prm
OPENMP-THREADS 12
digits 8
RESP-WEIGHT 1

#############################
##                         ##
##  Literature References  ##
##                         ##
#############################

Walker, B., Liu, C., Wait, E., Ren, P., J. Comput. Chem. 2022, 1. https://doi.org/10.1002/jcc.26954

Wu, J.C.; Chattree, G.; Ren, P.Y.; Automation of AMOEBA polarizable force field
parameterization for small molecules. Theor Chem Acc.

atom          401    401    C     "DMC                 "         6    12.011    3
atom          402    402    O     "DMC                 "         8    15.999    1
atom          403    403    O     "DMC                 "         8    15.999    2
atom          405    405    H     "DMC                 "         1     1.008    1
atom          404    404    C     "DMC                 "         6    12.011    4

# amoeba09 matching SMARTS from molecule  [['[#6](-[H])(-[H])-[H]', [1]]] to SMARTS from parameter file [['[#8](-[#6](-[H])(-[H])-[H])-[#6](-[H])(-[H])-[H]', [2]]] with tinker type descriptions [[('C', '"Methyl Ether CH3"')]]
# Missing vdw parameters
# [404] = [[9, 5]]
vdw 404 3.8200 0.1010

# amoeba09 matching SMARTS from molecule  [['[#6](-[H])(-[H])-[H]', [4]]] to SMARTS from parameter file [['[#8](-[#6](-[H])(-[H])-[H])-[#6](-[H])(-[H])-[H]', [5]]] with tinker type descriptions [[('H', '"Methyl Ether H3C"')]]
# Missing vdw parameters
# [405] = [[6, 7, 8, 10, 11, 12]]
vdw 405 2.8900 0.0240 0.910

# amoeba09 matching SMARTS from molecule  [['[#6](-[#8])=[#8]', [3]]] to SMARTS from parameter file [['[#8](-[#6](=[#8])-[#6](-[H])(-[H])-[H])-[H]', [3]]] with tinker type descriptions [[('O', '"Acetic Acid O=C"')]]
# Missing vdw parameters
# [402] = [[2]]
vdw 402 3.3000 0.1120

# amoeba09 matching SMARTS from molecule  [['[#6]-[#8]', [2]]] to SMARTS from parameter file [['[#8](-[#6](=[#8])-[#6](-[H])(-[H])-[H])-[H]', [1]]] with tinker type descriptions [[('O', '"Acetic Acid OH"')]]
# Missing vdw parameters
# [403] = [[3, 4]]
vdw 403 3.4050 0.1100

# amoeba09 matching SMARTS from molecule  [['[#6](-[#8])=[#8]', [1]]] to SMARTS from parameter file [['[#8](-[#6](=[#8])-[#6](-[H])(-[H])-[H])-[H]', [2]]] with tinker type descriptions [[('C', '"Acetic Acid C=O"')]]
# Missing vdw parameters
# [401] = [[1]]
vdw 401 3.7800 0.1060

# amoeba21 comments=O-est, Oxygen of carbonyl group (correct the O=C in ester) C=O, sp2 carbon, carbonic ester SMARTS match = [$([OH0]=[C]([O]([C])))] [CX3](=O)(O[C])(O[C])
# [403, 401] = [[3, 4], [1]]
bond 403 401 292.725366 1.34

# amoeba21 comments=O-est, Oxygen of carbonyl group (correct the O=C in ester) C=O, sp2 carbon, carbonic ester SMARTS match = [$([OH0]=[C]([O]([C])))] [CX3](=O)(O[C])(O[C])
# [402, 401] = [[2], [1]]
bond 402 401 292.725366 1.22

# amoeba21 comments=C33, sp3 carbon with 3 Hydrogens O-est, Oxygen of carbonyl group (correct the O=C in ester) SMARTS match = [CX4H3] [OX2H0]([C](=O))[C]
# [404, 403] = [[9, 5], [3, 4]]
bond 404 403 317.390148 1.44

# amoeba21 comments=HCH3, H on sp3 carbon, -CH3 C33, sp3 carbon with 3 Hydrogens SMARTS match = [H][CH3] [CX4H3]
# [405, 404] = [[6, 7, 8, 10, 11, 12], [9, 5]]
bond 405 404 345.97894 1.09

# amoeba21 comments=C33, sp3 carbon with 3 Hydrogens O-est, Oxygen of carbonyl group (correct the O=C in ester) C=O, sp2 carbon, carboxylic ester SMARTS match = [CX4H3] [OX2H0]([C](=O))[C] [$([CX3](=O)([O]([C])))]
# [404, 403, 401] = [[9, 5], [3, 4], [1]]
angle 404 403 401 65.992605 108.3131


# amoeba21 comments=O-est, Oxygen of carbonyl group (correct the O=C in ester) C=O, sp2 carbon, carbonic ester O-est, Oxygen of carbonyl group (correct the O=C in ester) SMARTS match = [OX2H0]([C](=O))[C] [CX3](=O)(O[C])(O[C]) [$([OH0]=[C]([O]([C])))]
# [403, 401, 402] = [[3, 4], [1], [2]]
angle 403 401 402 54.469344 126.24


# amoeba21 comments=O-est, Oxygen of carbonyl group (correct the O=C in ester) C=O, sp2 carbon, carbonic ester O-est, Oxygen of carbonyl group (correct the O=C in ester) SMARTS match = [OX2H0]([C](=O))[C] [CX3](=O)(O[C])(O[C]) [$([OH0]=[C]([O]([C])))]
# [403, 401, 403] = [[3, 4], [1], [3, 4]]
angle 403 401 403 54.469344 109.3521


# amoeba21 comments=HCH3, H on sp3 carbon, -CH3 C33, sp3 carbon with 3 Hydrogens O-est, Oxygen of carbonyl group (correct the O=C in ester) SMARTS match = [H][CH3] [CX4H3] [OX2H0]([C](=O))[C]
# [405, 404, 403] = [[6, 7, 8, 10, 11, 12], [9, 5], [3, 4]]
angle 405 404 403 60.030198 108.6


# amoeba21 comments=HCH3, H on sp3 carbon, -CH3 C33, sp3 carbon with 3 Hydrogens HCH3, H on sp3 carbon, -CH3 SMARTS match = [H][CH3] [CX4H3] [H][CH3]
# [405, 404, 405] = [[6, 7, 8, 10, 11, 12], [9, 5], [6, 7, 8, 10, 11, 12]]
angle 405 404 405 36.406852 110.32


# amoeba21 comments=C33, sp3 carbon with 3 Hydrogens O-est, Oxygen of carbonyl group (correct the O=C in ester) C=O, sp2 carbon, carboxylic ester SMARTS match = [CX4H3] [OX2H0]([C](=O))[C] [$([CX3](=O)([O]([C])))]
# [404, 403, 401] = [[9, 5], [3, 4], [1]]
strbnd 404 403 401 21.334 -21.334

# amoeba21 comments=O-est, Oxygen of carbonyl group (correct the O=C in ester) C=O, sp2 carbon, carbonic ester O-est, Oxygen of carbonyl group (correct the O=C in ester) SMARTS match = [OX2H0]([C](=O))[C] [CX3](=O)(O[C])(O[C]) [$([OH0]=[C]([O]([C])))]
# [403, 401, 402] = [[3, 4], [1], [2]]
strbnd 403 401 402 7.6289 7.6289

# amoeba21 comments=O-est, Oxygen of carbonyl group (correct the O=C in ester) C=O, sp2 carbon, carbonic ester O-est, Oxygen of carbonyl group (correct the O=C in ester) SMARTS match = [OX2H0]([C](=O))[C] [CX3](=O)(O[C])(O[C]) [$([OH0]=[C]([O]([C])))]
# [403, 401, 403] = [[3, 4], [1], [3, 4]]
strbnd 403 401 403 7.6289 7.6289

# amoeba21 comments=HCH3, H on sp3 carbon, -CH3 C33, sp3 carbon with 3 Hydrogens O-est, Oxygen of carbonyl group (correct the O=C in ester) SMARTS match = [H][CH3] [CX4H3] [OX2H0]([C](=O))[C]
# [405, 404, 403] = [[6, 7, 8, 10, 11, 12], [9, 5], [3, 4]]
strbnd 405 404 403 5.7126 5.7126

# amoeba21 comments=HCH3, H on sp3 carbon, -CH3 C33, sp3 carbon with 3 Hydrogens HCH3, H on sp3 carbon, -CH3 SMARTS match = [H][CH3] [CX4H3] [H][CH3]
# [405, 404, 405] = [[6, 7, 8, 10, 11, 12], [9, 5], [6, 7, 8, 10, 11, 12]]
strbnd 405 404 405 5.7126 5.7126

# amoeba21 comments=O-est, Oxygen of carbonyl group (correct the O=C in ester) C=O, sp2 carbon, carbonate SMARTS match = [OX2H0]([C](=O))[C] [CX3](=O)([OH0])[OH0]
# [403, 401] = [[3, 4], [1]]
opbend 403 401 0 0 36.3022

# amoeba21 comments=O-est, Oxygen of carbonyl group (correct the O=C in ester) C=O, sp2 carbon, carbonate SMARTS match = [OX2H0]([C](=O))[C] [CX3](=O)([OH0])[OH0]
# [402, 401] = [[2], [1]]
opbend 402 401 0 0 36.3022

# amoeba09 matching SMARTS from molecule  [['[*]~[*]~[*]~[*]', [1, 2, 3, 4]]] to SMARTS from parameter file [['[#6](-[H])(-[H])(-[H])-[#6](-[H])(-[H])(-[H])', [2, 1, 5, 6]]] with tinker type descriptions [[('H', '"Alkane H3C-"'), ('C', '"Alkane CH3-"'), ('C', '"Alkane CH3-"'), ('H', '"Alkane H3C-"')]]
# Missing torsion parameters, will attempt to fit parameters
# [405, 404, 403, 401] = [[6, 7, 8, 10, 11, 12], [9, 5], [3, 4], [1]]
# Fitted torsion
# Torsion 401 403 404 405 RMSD(MM2,QM) 0.029402340151312573 RelativeRMSD(MM2,QM) 0.019305616918499937 
 torsion     405  404  403  401      0.013 0.0 1   0.013 180.0 2   0.129 0.0 3

# amoeba09 matching SMARTS from molecule  [['[*]~[*]~[*]~[*]', [1, 2, 3, 4]]] to SMARTS from parameter file [['[#6](-[H])(-[H])(-[H])-[#6](-[H])(-[H])(-[H])', [2, 1, 5, 6]]] with tinker type descriptions [[('H', '"Alkane H3C-"'), ('C', '"Alkane CH3-"'), ('C', '"Alkane CH3-"'), ('H', '"Alkane H3C-"')]]
# Missing torsion parameters, will attempt to fit parameters
# [404, 403, 401, 402] = [[9, 5], [3, 4], [1], [2]]
# Fitted torsion
# Torsion 402 401 403 404 RMSD(MM2,QM) 0.12396100949787807 RelativeRMSD(MM2,QM) 0.014633198497850615 Boltzmann Fit
 torsion     404  403  401  402      0.357 0.0 1   3.418 180.0 2  -0.089 0.0 3

# amoeba09 matching SMARTS from molecule  [['[*]~[*]~[*]~[*]', [1, 2, 3, 4]]] to SMARTS from parameter file [['[#6](-[H])(-[H])(-[H])-[#6](-[H])(-[H])(-[H])', [2, 1, 5, 6]]] with tinker type descriptions [[('H', '"Alkane H3C-"'), ('C', '"Alkane CH3-"'), ('C', '"Alkane CH3-"'), ('H', '"Alkane H3C-"')]]
# Missing torsion parameters, will attempt to fit parameters
# [404, 403, 401, 403] = [[9, 5], [3, 4], [1], [3, 4]]
# Fitted torsion
# Torsion 403 401 403 404 RMSD(MM2,QM) 0.12396100949787807 RelativeRMSD(MM2,QM) 0.014633198497850615 Boltzmann Fit
 torsion     404  403  401  403     -0.543 0.0 1   3.418 180.0 2  -0.097 0.0 3

#SOLUTE-SMARTS 405 [#1][CH3]

SOLUTE 405 3.143 3.374 4.144

#SOLUTE-SMARTS 401 [#6;D3]([*])([*])(=[*])

SOLUTE 401 2.9301 4.506 4.2118

#SOLUTE-SMARTS 404 [CH3]

SOLUTE 404 3.5062 3.309 4.536

#SOLUTE-SMARTS 402 [#8;D1]=[#6]

SOLUTE 402 2.9835 3.356 2.9616

#SOLUTE-SMARTS 403 [O;D2][#6]

SOLUTE 403 3.1684 3.134 3.0999


# amoeba21 comments=C on carbonyl O SMARTS match = [CX3](=O)
# [401] = [[1]]
polarize           401          2.0730     0.3900 402 403
# amoeba21 comments=O on carbonyl group SMARTS match = [OX1]=[CX3]
# [402] = [[2]]
polarize           402          2.3740     0.3900 401
# amoeba21 comments=O connected to carbonyl group SMARTS match = [OX2][CX3]
# [403] = [[3, 4]]
polarize           403          2.3740     0.3900 401
# amoeba21 comments=H on non-aromatic carbon SMARTS match = [H][C]
# [405] = [[6, 7, 8, 10, 11, 12]]
polarize           405          0.4800     0.3900 404
# amoeba21 comments=C on non-aromatic oxygen SMARTS match = [C][O]
# [404] = [[9, 5]]
polarize           404          1.6200     0.3900 405

#
# Multipoles from Electrostatic Potential Fitting
#

# [401] = [[1]]
multipole   401 -403 -403               0.98274
                                        0.00000    0.00000    0.19848
                                        0.13103
                                        0.00000   -0.31919
                                        0.00000    0.00000    0.18816
# [402] = [[2]]
multipole   402  401  403              -0.66638
                                        0.00003    0.00000    0.02018
                                       -0.39219
                                        0.00000    0.15297
                                       -0.00035    0.00000    0.23922
# [403] = [[3, 4]]
multipole   403  401  404              -0.38755
                                        0.32060    0.00000    0.01030
                                        0.36329
                                        0.00000   -0.39779
                                       -0.20356    0.00000    0.03450
# [404] = [[9, 5]]
multipole   404  403  401               0.05936
                                       -0.04997    0.00000    0.30946
                                       -0.44459
                                        0.00000   -0.42084
                                       -0.03445    0.00000    0.86543
# [405] = [[6, 7, 8, 10, 11, 12]]
multipole   405  404  403               0.05667
                                       -0.01072    0.00000   -0.11717
                                        0.02360
                                        0.00000    0.03369
                                       -0.02617    0.00000   -0.05729

代码

文本

poltype2还对目标分子的一些关键参数，主要是二面角进行了拟合结果和QM方法的比较，可以在OPENME文件夹中找到图片

代码

文本

[15]

from IPython.display import Image

for image in Path('OPENME').glob('*png'):

print(image.stem)

display(Image(image))

print(Path('OPENME/README.txt').read_text())

DMC-energy-1-3

DMC-energy-1-3_use_weights

DMC-energy-3-9

DMC-fit-1-3

DMC-fit-1-3_useweights

DMC-fit-3-9

Please ensure that the Root Mean Square Deviation (RMSD) between QM (red curve) and MM2 post-fit torsion (blue curve) has a decent fit.
Plots with *energy* contain the QM energy vs dihedral angle vs MM2 energy vs dihedral angle
Plots with *fit* contain the QM-MM1 (prefit torsion energy) vs dihedral angle and the fit spline
The first two numbers in plot are the rotatable bond in the parent, last two are the rotatable bond in the fragment
Plots with _use_weights in the filename are Boltzman fitted plots
Please ensure that the QM dipole and MM dipole match reasonably well
Please ensure that the RMSD between average QM potential and average MM potential is small
Tue Aug 22 00:54:03 2023 Absolute or Relative RMSPD of QM and MM torsion profiles is high, RMSPD = 1.0747848705527308 Tolerance is 1 kcal/mol RMSPDRel =0.12813249099302118 tolerance is 0.1
Tue Aug 22 00:56:10 2023 Relative error of 0.012396694214876044 for QMDipole 0.242 and 0.239 for MMDipole  tolerance = 0.5 /home/input_lbg-13211-8456961/Temp

代码

文本

补充

该文件可以直接用于Tinker进行MD（但要注意修改原子类型的定义），也可以借助OpenMM使用GPU加速运行AMOEBA力场，OpenMM 需要专门的xml格式力场文件，您可以使用随数据集一并提供的 read_poltype2_output.py文件转化获得xml文件。这里不再赘述，可以参考http://docs.openmm.org/latest/userguide/library/08_amoeba_plugin.html 。

代码

文本

参考资料

https://github.com/TinkerTools/poltype2
Yue Shi, Zhen Xia, Jiajing Zhang, Robert Best, Chuanjie Wu, Jay W. Ponder, and Pengyu Ren. Journal of Chemical Theory and Computation 2013 9 (9), 4046-4063 DOI: 10.1021/ct4003702 [ https://doi.org/10.1021/ct4003702 ]
Jay W. Ponder, Chuanjie Wu, Pengyu Ren, Vijay S. Pande, John D. Chodera, Michael J. Schnieders, Imran Haque, David L. Mobley, Daniel S. Lambrecht, Robert A. DiStasio Jr., Martin Head-Gordon, Gary N. I. Clark, Margaret E. Johnson, and Teresa Head-Gordon. The Journal of Physical Chemistry B 2010 114 (8), 2549-2564 DOI: 10.1021/jp910674d [ https://doi.org/10.1021/jp910674d ]

代码

文本

分子动力学

力场

拟合

分子动力学力场拟合

已赞1

本文被以下合集收录

镜影映照

更新于 2024-04-29

1 篇0 人关注