sudo yum install epel-releasesudo yum install -y python-ase安装ASE之前必须安装python和numpy，如果想用图形界面或是画图处理，那么可以安装matplotlib，pygtk等。
具体要求可以查看ASE installnation requirements。
C:\Users\Tang\AppData\Local\Enthought\Canopy\User\Scripts添加到系统环境变量Path里。Win+R打开cmd，输入pip install ase，这样ase就安装好了，你可以用它来处理一些数据，但无法使用ase-gui。如果想运行ase-gui，那就必须安装win-64位的GTK+和pygtk。
C:\opt\gtk，然后把C:\opt\gtk\bin添加到系统环境变量Path里，你可以在cmd里输入gtk-demo来测试。cmd，输入pip install pycairo_gtk-1.10.0-cp27-none-win_amd64.whl来安装，剩下两个文件同理。cmd中输入ase-gui，可以看到ASE的图形界面。pygtk的安装如果有问题，可以参考python pygtk windows 7 64 bit - stackOverflow。 如果安装没有成功，只是不能用图形界面而已，不影响ASE里面的一些函数的使用。
结构参数可以从免费的晶体数据库找到American Mineralogist Crystal Structure Database
在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])
切一个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)
金属的表面可以直接切得，比如Pd(111)面可以这样得到
from ase.lattice.surface import fcc111 slab = fcc111('Pd', a=3.977, size=(2,2,3), vacuum=10.0, orthogonal=True)
类似的方法还有fcc100，fc110，bcc100，hcp0001，hcp10m10，diamond100，diamond111等。
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])
构建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)
添加吸附物有两种方法，一种是用内置的add_adsorbate函数，一种是手动添加。
把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)
把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")
atoms.center(10, axis=2)：先居中，然后在Z方向上，上下对称地加10A真空层，一共是20A的真空层。以建一个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")
atoms.set_constraint()atoms = atoms * (a, b, c)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/Å
MinModeTranslate(Optimizer) has logged twice, one in MinModeTranslate, another in Optimizer.from ase.cluster import Icosahedron from ase.calculators.lj import LennardJones from ase.md.velocitydistribution import MaxwellBoltzmannDistribution from ase.md.nvtberendsen import NVTBerendsen from ase.units import kB, fs temperature = 1000 atoms = Icosahedron('C', 2, 1.53) calc = LennardJones() atoms.calc = calc MaxwellBoltzmannDistribution(atoms, temperature_K=temperature) dyn = NVTBerendsen(atoms, 1*fs, temperature_K=temperature, taut=100*fs, logfile='md.log', trajectory='md.traj', loginterval=100 ) # 1 ns dyn.run(1e6)
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)
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__().
测试脚本，提交任务的命令为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())
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
[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