Atomistic Simulation Environment
Table of Contents

Resources

Compatibility

Installation

Centos

Windows(64-bit)

安装ASE之前必须安装pythonnumpy,如果想用图形界面或是画图处理,那么可以安装matplotlibpygtk等。 具体要求可以查看ASE installnation requirements

scientific python

ase-gui

如果想运行ase-gui,那就必须安装win-64位的GTK+和pygtk。

pygtk的安装如果有问题,可以参考python pygtk windows 7 64 bit - stackOverflow。 如果安装没有成功,只是不能用图形界面而已,不影响ASE里面的一些函数的使用。

Building model

Bulk

结构参数可以从免费的晶体数据库找到American Mineralogist Crystal Structure Database

TiO2

American Mineralogist Crystal Structure Database搜索anatase得到参数,a和c是优化好的晶格参数。 如何优化可以参考用VASP优化TiO2_anatase晶格常数(a, c)

from ase.lattice.spacegroup import crystal

a = 3.862
c = 9.551
TiO2 = crystal(['Ti', 'O'], basis=[(0 ,0, 0), (0, 0, 0.208)], spacegroup=141, cellpar=[a, a, c, 90, 90, 90])

Surface

common surface

切一个TiO2的101面,层数为4,上下均加5A的真空层。对于TiO2的101面,不同的终端稳定性不同。 最稳定的终端可以从文献里找到,比如这篇Structure and energetics of stoichiometric TiO2 anatase surfaces。 如果你想得到4层指定终端的101面,直接用下面的命令是不行的。一般需要多切几层,比如6层或8层,然后把多余的原子删掉,就可以得到指定终端的4层的101面

from ase.lattice.surface import surface

slab101 = surface(TiO2, (1, 0, 1), 4, vacuum=5)

metal surface

金属的表面可以直接切得,比如Pd(111)面可以这样得到

from ase.lattice.surface import fcc111
slab = fcc111('Pd', a=3.977, size=(2,2,3), vacuum=10.0, orthogonal=True)

类似的方法还有fcc100fc110bcc100hcp0001hcp10m10diamond100diamond111等。

sqrt surface

from ase.build import fcc111, cut
from ase.visualize import view

slab = fcc111('Cu', (1, 1, 4), vacuum=10)

slab_1_sqrt3 = cut(slab, (1, 0, 0), (-1, 2, 0), (0, 0, 1))
# a, b, c: cell vectors
# new_a = old_a
# new_b = -old_a + 2*old_b
# new_c = old_c

slab_sqrt3_1 = cut(slab, (2, -1, 0), (0, 1, 0), (0, 0, 1))
slab_sqrt3_sqrt3 = cut(slab, (1, 1, 0), (-1, 2, 0), (0, 0, 1))
view([slab, slab_1_sqrt3, slab_sqrt3_1, slab_sqrt3_sqrt3])

2-D materials with hexagonal structure

构建3层六角形的二维材料,例如VS2。支持2H和1T两种类型

from ase.lattice.surface import mx2
vs2 = mx2(formula='VS2', kind='1T', a=3.20, thickness=2.852, size=(1, 1, 1), vacuum=10)

Atoms operations

add adsorbate

添加吸附物有两种方法,一种是用内置的add_adsorbate函数,一种是手动添加。

add_adsorbate function

把O2加到Pd(100)的bridge位上可以像下面的例子这样做,先搭建O2和Pd(100)的模型,然后找到吸附位置的坐标。(坐标需要view(Pd100)查看) 这里需要注意的是,ase在添加的时候是把O2的第一个坐标放到指定的吸附坐标上,所以想把O2的中心放到bridge的中心,吸附的坐标应为(bridge_middle[0]-r, bridge_middle[1])。 吸附完之后,再把整体放到正中央,并对称地加上真空层。

from ase.all import *
from ase.lattice.surface import fcc100, add_adsorbate
from ase.visualize import view

d = 1.24414
r = d/2 # O2 radius
h = 3 # adsorbate height

# Build Pd100 and O2
Pd100 = fcc100("Pd", a=3.977, size=(4,4,5), vacuum=10)
O_2 = Atoms("2O", positions=[(0, 0, 0), (d, 0, 0)])

# Get adsorbate position
bridge_Pd = [Pd100.get_positions()[73], Pd100.get_positions()[74]]
bridge_middle = ((bridge_Pd[0][0]+bridge_Pd[1][0])/2, (bridge_Pd[0][1]+bridge_Pd[1][1])/2)

# add O2 on Pd100
add_adsorbate(Pd100, O_2, h, (bridge_middle[0]-r, bridge_middle[1]))

# Make the whole model center, and add 10 vaccum on Z-direction
Pd100.center(10, axis=2)

# Check the model with ase-gui
view(Pd100)

add adsorbate manually

把Pd(111)加到TiO2(101)上面,TiO2(101)从优化好的CONTCAR中读取。加的时候先把Pd(111)的晶格设置成TiO2(101)一样,然后放到中间。 这个时候,使用+=,Pd(111)会跑到TiO2(101)里面,所以需要把Pd(111)。调高的高度根据Pd最下面的z坐标和TiO2最上面的z坐标而定。

from ase.all import *
from ase.utils.geometry import *
from ase.lattice.surface import *
from ase.visualize import view

Pd111 = fcc111("Pd", a=3.977, size=(3,4,3), orthogonal=False)

slab101 = read("CONTCAR", format="vasp")


# Add Pd111 to slab101 
cell = slab101.cell
Pd111.set_cell(cell)

Pd111.center()

for atom in Pd111:
  atom.z += 3.34

slab101 += Pd111

#adsorb_site = [slab101.cell[0][0], slab101.cell[1][1]]
#add_adsorbate(slab101, Pd111, 2, adsorb_site, mol_index=2)

# center slab101 and add vaccum
slab101.center(10, 2)

# fix atoms
fixed_atoms = FixAtoms(indices=[atom.index for atom in slab101 if atom.z < 14])
slab101.set_constraint(fixed_atoms)

# sort atoms
slab101 = sort(slab101)

view(slab101)

#slab101.write("POSCAR")

center atoms

delete atoms

以建一个4层TiO2_anatase为例,来介绍ASE中删除原子的做法。根据原子的index选择需要删除的原子,直接del

import ase
from ase.all import *
from ase.utils.geometry import *
from ase.lattice.spacegroup import crystal
from ase.lattice.surface import surface, add_vacuum
from ase.visualize import view

a = 3.862
c = 9.551
TiO2 = crystal(['Ti', 'O'], basis=[(0 ,0, 0), (0, 0, 0.208)], spacegroup=141, cellpar=[a, a, c, 90, 90, 90])

slab101 = surface(TiO2, (1, 0, 1), 6, vacuum=6)

# 1x3x1 supercell
slab101 = slab101*(1, 3, 1)

# delete atoms
extra_atoms = []
for atom in slab101:
    if atom.z > 22 or atom.z < 8:
        extra_atoms.append(atom.index)
del slab101[extra_atoms]

# fix atoms
fixed_atoms = FixAtoms(indices=[atom.index for atom in slab101 if atom.z < 15])
slab101.set_constraint(fixed_atoms)

# sort atoms
slab101 = sort(slab101)

view(slab101)
slab101.write("POSCAR")

fix atoms

supercell

NEB

import numpy as np
from ase.io import read
from ase.neb import NEB

trajs = read('neb_rough.traj', index=':7')

neb = NEB(images=trajs, method='eb')
neb_f = neb.get_forces()
neb_max_f = np.sqrt((neb_f**2).sum(axis=1).max())
neb_e = neb.get_potential_energy()

print(f'neb        {neb_e:.4f} eV      {neb_max_f:.4f} eV/Å')

for i in range(len(trajs)):
    atoms = trajs[i]
    e = atoms.get_potential_energy()
    f = atoms.get_forces()
    f_max = np.sqrt((f**2).sum(axis=1).max())
    print(f'image-{i}    {e:.4f} eV      {f_max:.4f} ev/Å')

# neb        -2712.6466 eV      6.5730 eV/Å     !!!!
#
# image-0    -2728.9443 eV      0.0187 ev/Å
# image-1    -2721.9550 eV      6.7862 ev/Å
# image-2    -2715.9325 eV      6.5309 ev/Å
# image-3    -2714.0065 eV      3.9710 ev/Å     
# image-4    -2712.6466 eV      4.0747 ev/Å     !!!!
# image-5    -2722.8021 eV      6.4575 ev/Å
# image-6    -2731.2418 eV      0.0168 ev/Å

Dimer

Plot

colors

from ase.data.colors import jmol_colors
from ase.visualize.plot import plot_atoms

from matplotlib.colors import to_rgb
colors = jmol_colors[atoms.get_atomic_numbers()]
colors[[37, 43, 45, 47]] = to_rgb('C1')  #

ax = plot_atoms(atoms, ax, colors=colors)

Calculator

CP2K

The first step is to setup cp2k environment in bash script. Remember to use cp2k_shell.xxx which is an interactive shell.

export ASE_CP2K_COMMAND="mpirun -n 16 /home/tang/opt/cp2k/cp2k-7.1.0/exe/Linux-x86-64-intel-minimal/cp2k_shell.popt"

The interactive shell behaves somewhat different than normal binary. It will keep alive if you don't kill the shell by hand. This could be a problem if you want to perform single point calculations on multiple structures. For example, the loop will not continue after the first calculation is finished because the CP2K shell is not terminated at the end of loop.

# not working
for atoms in traj:
  calc = CP2K()
  atoms.set_calculator(calc)
  atoms.get_forces()

The easiest solution is to define the calculator class before the loop so that all steps within the loop will use the same CP2K shell.

# not working
calc = CP2K()
for atoms in traj:
  atoms.set_calculator(calc)
  atoms.get_forces()

The problem here is that energy and force information are erased when the same calculator is attached to different atoms. On has to terminate CP2K shell by hand either via calc.__del__() or atoms.calc.__del__().

# working
for atoms in traj:
  calc = CP2K()
  atoms.set_calculator(calc)
  atoms.get_forces()

  calc.__del__()

del calc or del atoms.calc will not call __del__().

ASE scripts

FUNCOM

test script

测试脚本,提交任务的命令为qsub vasp.py

#!/home/tangzeyuan/bin/python2.7
#$ -S /home/tangzeyuan/bin/python2.7
#$ -cwd
#$ -j y
#$ -N ase_test
#$ -pe make 4
import sys
import os
from subprocess import Popen, PIPE
from ase import Atoms, Atom
from ase.calculators.vasp import Vasp

os.environ["VASP_COMMAND"] = 'mpirun -r ssh -np 4 /home/tangzeyuan/bin/vasp'
os.environ["VASP_PP_PATH"] = '/home/tangzeyuan/lib/vasp'

### Prepare MPI Environment (Do not modify codes below !!!) ###
### apply to 42.244.24.64                   ###
def source(script, update=1):
    pipe = Popen(". %s; env" % script, stdout=PIPE, shell=True)
    data = pipe.communicate()[0]

    env = dict((line.split("=", 1) for line in data.splitlines()))
    if update:
        os.environ.update(env)

    return env
source("/share/apps/intel/Compiler/11.1/073/bin/iccvars.sh intel64")
source("/share/apps/intel/Compiler/11.1/073/bin/ifortvars.sh intel64")
source("/share/apps/intel/impi/3.2.0.011/bin64/mpivars.sh")
### Prepare MPI Environment (Do not modify codes above !!!) ###

### Build model
a = [6.5, 6.5, 7.7]
d = 2.3608
NaCl = Atoms([Atom('Na', [0, 0, 0], magmom=1.928),
              Atom('Cl', [0, 0, d], magmom=0.75)],
              cell=a)

### Set up INCAR, KPOINTS, POTCAR
calc = Vasp(prec='Accurate',
            xc='PBE',
            lreal=False,
            lcharg=False,
            gamma=True, # use gama-sampling instead of Monkhorst-Pack
            kpts=(1, 1, 1))

### Set up input files & run vasp
NaCl.set_calculator(calc)

### Data processing after vasp finishing
print(NaCl.get_magnetic_moment())

encut test

Pd的encut测试

#!/home/tangzeyuan/bin/python2.7
#$ -S /home/tangzeyuan/bin/python2.7
#$ -cwd
#$ -j y
#$ -N ase_test
#$ -pe make 12
import sys
import os
from subprocess import Popen, PIPE
from ase import Atoms, Atom
from ase.lattice import bulk
from ase.calculators.vasp import Vasp

os.environ["VASP_COMMAND"] = 'mpirun -r ssh -np 12 /home/tangzeyuan/bin/vasp'
os.environ["VASP_PP_PATH"] = '/home/tangzeyuan/lib/vasp'

### Prepare MPI Environment (Do not modify codes below !!!) ###
### apply to 42.244.24.64                   ###
def source(script, update=1):
    pipe = Popen(". %s; env" % script, stdout=PIPE, shell=True)
    data = pipe.communicate()[0]

    env = dict((line.split("=", 1) for line in data.splitlines()))
    if update:
        os.environ.update(env)

    return env
source("/share/apps/intel/Compiler/11.1/073/bin/iccvars.sh intel64")
source("/share/apps/intel/Compiler/11.1/073/bin/ifortvars.sh intel64")
source("/share/apps/intel/impi/3.2.0.011/bin64/mpivars.sh")
### Prepare MPI Environment (Do not modify codes above !!!) ###

Pd = bulk('Pd', 'fcc', a=3.89)

encut_list = [200, 300, 400, 500, 600]
energy_list = []

for en in encut_list:
    calc = Vasp(prec='Accurate',
                xc='PBE',
                lreal=False,
                lcharg=False,
                encut=en,
                gamma=True, # use gama-sampling instead of Monkhorst-Pack
                kpts=(8, 8, 8))
    Pd.set_calculator(calc)
    energy_list.append(Pd.get_potential_energy())

print encut_list
print energy_list

FAQ

[tangzeyuan@hn 01.MoSe2]$ pip install --upgrade ase
Collecting ase
  Could not fetch URL https://pypi.python.org/simple/ase/: There was a problem confirming the ssl certificate: [SSL: CERTIFICATE_VERIFY_FAILED] certificate verify failed (_ssl.c:590) - skipping
  Could not find a version that satisfies the requirement ase (from versions: )
No matching distribution found for ase

解决办法

[tangzeyuan@hn 01.MoSe2]$ pip --trusted-host pypi.python.org install --upgrade ase 
Collecting ase
  Downloading ase-3.10.0.tar.gz (1.6MB)
    100% |################################| 1.6MB 10.5MB/s 
Installing collected packages: ase
  Running setup.py install for ase ... done
Successfully installed ase-3.10.0