Pysimm教程——从Examples解读和学习建模和计算01_methane甲烷

我搜寻过B站还有各种渠道，甚至是Pysimm官网，都没有一个系统性的教程关于pysimm。

我今天就来抛出这块砖，请各位大佬批评指正：

所出现的代码版权和解释权都归于：

pysimm该网站所有

 且先看create.py这个文件和README.md这两个原始文件

create.py:

from pysimm import system, lmps, forcefield

def run(test=False):
    # 创建一个空系统
    print('Example progress: Creating an empty system...')
    s = system.System()
    
    # 在我们的系统中创建一个新分子
    print('Example progress: Adding an empty molecule container to our system...')
    m = s.molecules.add(system.Molecule())
    
    # 检索 GAFF2 参数
    print('Example progress: Retrieving Dreiding force field parameters...')
    f = forcefield.Gaff2()
    s.forcefield = f.name
    
    # get a copy of the c3 particle type object from GAFF 找到一个C3粒子类型
    # get method returns a list, we need the first element 需要第一个元素
    gaff_c3 = s.particle_types.add(f.particle_types.get('c3')[0].copy())
    
    # get hc particle type object from GAFF 得到hc粒子类型从GAFF中
    gaff_hc = s.particle_types.add(f.particle_types.get('hc')[0].copy())
    
    # we'll first make the carbon atom at the origin 原始中将C作为第一个C原子
    # we'll include gasteiger charges later  后续会给予电荷
    print('Example progress: Adding carbon atom at origin...')
    # 将第一个C原子放在了原点，电荷为0
    c1 = s.particles.add(system.Particle(type=gaff_c3, x=0, y=0, z=0, charge=0, molecule=m))
    
    # now we'll add 4 hydrogen atoms bonded to our carbon atom
    # 现在给C原子添加四个H原子
    # these atoms will be placed randomly 1.5 angstroms from the carbon atom
    # 这些原子随机的以1.5埃距离C的位置摆放
    # we'll optimize the structure using LAMMPS afterwords
    # 随后将使用LAMMPS进行结构优化
    # we supply the GAFF forcefield object so that bond and angle types can be added as well
    # 随后施加GAFF力场，添加bond和angle类型
    print('Example progress: Adding 4 hydrogen atoms at random positions bonded to the carbon atom...')
    h1 = s.add_particle_bonded_to(system.Particle(type=gaff_hc, charge=0, molecule=m), c1, f)
    h2 = s.add_particle_bonded_to(system.Particle(type=gaff_hc, charge=0, molecule=m), c1, f)
    h3 = s.add_particle_bonded_to(system.Particle(type=gaff_hc, charge=0, molecule=m), c1, f)
    h4 = s.add_particle_bonded_to(system.Particle(type=gaff_hc, charge=0, molecule=m), c1, f)
    
    # let's add gasteiger charges
    # 添加gasteiger电荷
    print('Example progress: Deriving Gasteiger charges...')
    s.apply_charges(f, charges='gasteiger')
    
    # right now there is no simulation box defined
    # 目前没有模拟box被定义
    # we'll define a box surrounding our methane molecule with a 10 angstrom padding
    # 我们定义了一个盒子，包围甲烷分子，10埃大小
    print('Example progress: Constructing Simulation box surrounding our new molecule...')
    s.set_box(padding=10)
    
    # before we optimize our structure, LAMMPS needs to know what type of 
    # pair, bond, and angle interactions we are using
    # these are determined by the forcefield being used
    # 在我们优化结构时，LAMMPS需要知道使用的相互作用，这里我们需要pair, bond和angle相互作用
    # 这些相互作用由使用的力场决定
    s.pair_style='lj'
    s.bond_style='harmonic'
    s.angle_style='harmonic'
    
    # we'll perform energy minimization using the fire algorithm in LAMMPS
    # 使用LAMMPS中的fire算法进行能量最小化
    print('Example progress: Optimizing structure using LAMMPS...')
    lmps.quick_min(s, min_style='fire', name='fire_min', etol=1e-10, ftol=1e-10)
    
    # write xyz, YAML, LAMMPS data, and chemdoodle json files
    print('Example progress: Saving structure to files...')
    s.write_xyz('methane.xyz')
    s.write_yaml('methane.yaml')
    s.write_lammps('methane.lmps')
    s.write_chemdoodle_json('methane.json')
    
    print('Example progress: Complete!')
    
if __name__ == '__main__':
    run()

其实真正看完python建模后，可以理解到其建模的思想，和鲍老师教的moltemplate建模思想一样，本教程的好处就是直接可以通过GAFF将模型计算出来。

下面我们好好研究一下人家的READme主要讲了什么吧：

### 引入 pysimm 模块/包

from pysimm import system, lmps, forcefield

####warning
###当你遇到**"ImportError: No module named pysimm"**确保你的pysimm路径在你的PYTHONPATH中

###准备一个空系统
##首先我们创建一个空的**system.System**项目并且将这个变量存储在**s**中。
##此系统项会包含并且管理所有的分子数据：

s = system.System()

##接下来我们需要在系统中添加分子。
##默认情况下，我们的系统 **s** 拥有一个包含项 **s.molecules**，我们需要将我们的新分子加入到里面
##我们创建一个新分子项目**system.Molecule()**并且传递给分子包含分类方法##**s.molecules.add()**
##此函数将新添加的项目返回至包含中，并且我们存储此变量为**m**

m = s.molecules.add(system.Molecule())

###现在我们拥有了一个包含所有分子数据的地方，我们需要获得描述力场原子相互作用的数据。
##一个**forcefield.Forcefield**项包含了计算你的系统的能量的必要参数以及逻辑性的输入规则
##从而分配原子类型于你的系统中。这个case中，我们将使用Gaff2 force field并且示例一个Gaff2力场项##并且将其存储在变量**f**中：

f = forcefield.Gaff2()

##此例假设使用者已经知道了Gaff2的原子类型，并且知道甲烷需要什么原子类型：c3和hc
##从**system.ParticleType**项中获得表示原子类型的**forcefield.Gaff2**

##forcefield.Gaff2中的项 f 包含一个项为 f.particle_types,这一项组织 ##system.ParticleType 项并且提供了解法 f.particle_types.get() 来提取特定的类型基于name
##此函数回到了一个列表 system.ParticleType项，但在此例中我知道只有一个项目被返回，所以我们得到##了list中的第一个元素。
## 为了取代添加 system.ParticleType 从 f 中，并且传递这个新创建的项目 system.ParticleType
## 至一个方法并且将其添加到粒子类型包含项目 s.particle_types。
##s.particle_types.add() 方法返回新添加的项，这个项我们作为Gaff2原子类型项储存

gaff_c3 = s.particle_types.add(f.particle_types.get('c3')[0].copy())
gaff_hc = s.particle_types.add(f.particle_types.get('hc')[0].copy())


###  添加原子进入system/molecule
##   目前我们有个系统项，s，一个分子项存储在我们的，m，两个Gaff2原子类型项，gaff_c3和gaff_hc
##   现在开始添加原子

### 首先，我们初始创建一个C原子。我们目前添加的粒子是0电荷的，后续会使用Gasteiger电荷。
## 示例一个 system.Particle 项目使用关键词 x, y, z, q，which part, particle tpye,
## gaff_c3

c1 = s.particles.add(system.Particle(type=gaff_c3, x=0, y=0, z=0, charge=0, molecule=m))

## sysgtem.System 拥有一个将新粒子连接与已经存在的粒子上的方法，前提是系统中已经有的粒子
## system.System.add_particle_bonded_to

h1 = s.add_particle_bonded_to(system.Particle(type=gaff_hc, charge=0, molecule=m),c1, f)

h2 = s.add_particle_bonded_to(system.Particle(type=gaff_hc, charge=0, molecule=m),c1, f)

h3 = s.add_particle_bonded_to(system.Particle(type=gaff_hc, charge=0, molecule=m),c1, f)

h4 = s.add_particle_bonded_to(system.Particle(type=gaff_hc, charge=0, molecule=m),c1, f)

## 后面我们可以将Gasteiger电荷派生使用 system.System.apply_charges() 类别办法
## 提供了forcefield.Forcefield 项并且证明我们希望使用Gasteiger approach

s.apply_charges(f, charge='gasteiger')

### 设置一个系统box并且优化结构
## 我们已经给系统添加了粒子，但是没有box. system.System 类别拥有能够建造粒子周围的box的能力
## 可选择性的，我们添加了10埃的box

s.set_box(padding=10)

## 我们优化结构的时候，LAMMPS需要知道我们使用的pair， bond， angle类型。归到system.System项

s.pair_style='lj'
s.bond_style='harmonic'
s.angle_style='harmonic'

## 我们使用两级minimize过程，首先使用最快的decent算法共轭梯度法
## lmps 模块在pysimm中包含快速的模拟方法
## 此case我们使用 lmps.quick_min()来传递我们的系统s和min_style我们想要用的方法
## 并且给每个模拟一个name以至于log文件中有可以识别的names

lmps.quick_min(s, min_style='sd', name='min_sd')
lmps.quick_min(s, min_style='cg', name='min_cg')

### 将系统写出为各种文件格式

## system.System 拥有各种系统data形式写成不同的file type。我们写出xyz, yaml, lammps data
## 和chemdoodle json 文件形式

s.write_xyz('methane.xyz')
s.write_yaml('methane.yaml')
s.write_lammps('methane.lmps')
s.write_chemdoodle_json('methane.json')