免费第一性原理软件的溶剂模型测试笔记

ultramanm87 · 发表于 Post on yesterday 08:55

本帖最后由 ultramanm87 于 2024-11-23 01:51 编辑

目的：作为一个搞环境的（重金属方向），如果要按照本论坛大佬们的建议，应该对体系（比如clay+金属离子）这种体系同时使用显式的、尤其是阳离子周围1-2层的水分子以及隐式溶剂模型，还应该进行一定的结构优化。之前搞cp2k sccs搞得头皮发麻尤其是想用它来优化结构，于是想试下其他软件是否能够在隐式溶剂模型下高效的结构优化，在这里做个笔记。欢迎讨论！

一、包含软件
CP2K
主要为高斯平面波基组，使用了SCCS模型(Andreussi et al., 2012)。实际计算中很依赖调整SCCS模型中判断溶质边界的参数，或者（据说）提高截断能，不然很容易不收敛或者能量振荡严重。
Quantum ESPRESSO
平面波软件，使用了SCCS模型(Andreussi et al., 2012)，需要使用外部的Environ插件来完成计算。默认参数够用。但是PW方法在原子数较多时计算实在太慢了。
GPAW
有PW、FD和LCAO三种方法。由于后续笔者想要在比较大的体系中进行计算，因此不考虑计算速度较慢的PW和FD方法。使用了CSM模型(Held and Walter, 2014)。
ABACUS
有PW和LCAO两种方法，同样PW方法在大体系很昂贵，因此只测试LCAO方法。使用了一种应用到VASP中的溶剂化方法(不是VASPSOL)(Mathew et al., 2014)。

二、小分子测试
测试平台7950x+64GB（no GPU）
将乙醇分子在隐式水溶剂中进行优化，比较优化的步数、优化的结构（pyscf的在非周期条件下，用B3LYP-D3(BJ)/def2-TZVP+IEF-PCM进行优化，带density fitting）、和优化完的溶解自由能（为省事儿使用溶剂模型中优化后的结构计算气象能量）。
能量参考:

以下是软件安装和测试的步骤。

GPAW 安装

GAPW可以直接用conda安装，安装命令为：

conda install conda-forge::gpaw

复制代码

上述命令建议先创建一个conda环境后再安装，避免环境冲突。使用上述安装后，可使用gpaw info查看有哪些安装的包：

| python-3.12.7 /home/haha/miniconda/envs/gpaw/bin/python |
| gpaw-24.6.0 /home/haha/miniconda/envs/gpaw/lib/python3.12/site-packages/gpaw/ |
| ase-3.23.0 /home/haha/miniconda/envs/gpaw/lib/python3.12/site-packages/ase/ |
| numpy-1.26.4 /home/haha/miniconda/envs/gpaw/lib/python3.12/site-packages/numpy/ |
| scipy-1.14.1 /home/haha/miniconda/envs/gpaw/lib/python3.12/site-packages/scipy/ |
| libxc-6.2.2 yes |
| _gpaw /home/haha/miniconda/envs/gpaw/lib/python3.12/site-packages/_gpaw.cpython-312-x86_64-linux-gnu.so |
| MPI enabled yes |
| OpenMP enabled yes |
| GPU enabled no |
| GPU-aware MPI no |
| CUPY /home/haha/miniconda/envs/gpaw/lib/python3.12/site-packages/gpaw/gpu/cpupy/__init__.py |
| scalapack yes |
| Elpa yes; version: 20211125 |
| FFTW yes |
| libvdwxc yes |
| PAW-datasets (1) /home/haha/miniconda/envs/gpaw/share/gpaw

复制代码

可见除了GPU相关，其他需要的库都已经装好了，可以直接使用mpirun -n 4 gpaw python script.py运行

模型构建
GPAW与ASE十分紧密，因此可以使用ASE创建分子与晶胞。创建ethanol_relax.py:

from ase.build import molecule
from ase.io import read as ase_read
from ase.io import write as ase_write
from ase.units import mol, kJ, kcal, Pascal, m
from ase.data.vdw import vdw_radi
from gpaw.convergence_criteria import Energy, Density, Eigenstates
from dftd3.ase import DFTD3 #需要安装simple-dftd3包
from ase.calculators.mixing import SumCalculator
from ase.optimize import BFGSLineSearch
from gpaw.solvation import (
# calculator
SolvationGPAW,
# cavities
EffectivePotentialCavity,
ADM12SmoothStepCavity,
#FG02SmoothStepCavity,
# custom classes for the cavities
Power12Potential,
ElDensity,
#SSS09Density,
# dielectric
LinearDielectric,
# non-el interactions
SurfaceInteraction,
VolumeInteraction,
LeakedDensityInteraction,
# surface and volume calculators
GradientSurface,
KB51Volume,
)#这边把所有主要的溶剂计算模块载入了，使用默认的CSM不需要这么多
atoms = molecule('CH3CH2OH')
atoms.center(vacuum=5)
atoms.pbc = [True, True, True]
ase_write("ethanol_initial.cif", atoms,format='cif') #保存初始结构给后续其他软件使用
h = 0.1 #格点精度
# solvent parameters for water from J. Chem. Phys. 141, 174108 (2014)
u0 = 0.180 # eV
epsinf = 78.36 # dimensionless
gamma = 18.4 * 1e-3 * Pascal * m # convert from dyne / cm to eV / Angstrom**2
T = 298.15 # Kelvin
p = -0.35 * 1e9 * Pascal
vdw_radii = vdw_radii.copy()
vdw_radii[1] = 1.09
vdw_radii[25] = 2.04
kappa_T = 4.53e-10 / Pascal
#需要用一个函数为原子添加半径参数
def atomic_radii(atoms_vdwr):
return [vdw_radii[n] for n in atoms_vdwr.numbers]
dft = SolvationGPAW(mode='lcao', basis='dzp', xc='PBE', nbands='110%', kpts=(2, 2, 2),
h=h, txt='ethanol_relax.txt', #输出文件
parallel=dict(band=2, # band parallelization
augment_grids=True, # use all cores for XC/Poisson
sl_auto=True, # enable ScaLAPACK parallelization
use_elpa=True),# enable elpa parallelization
cavity=EffectivePotentialCavity(
effective_potential=Power12Potential(atomic_radii, u0),
temperature=T,
surface_calculator=GradientSurface(),
volume_calculator=KB51Volume(compressibility=kappa_T, temperature=T)),
dielectric=LinearDielectric(epsinf=epsinf),
interactions=[SurfaceInteraction(surface_tension=gamma),
VolumeInteraction(pressure=p),
LeakedDensityInteraction(voltage=1)])
atoms.calc = SumCalculator([DFTD3(method="PBE", damping="d3bj"), dft]) #加上色散校正
relax = BFGSLineSearch(atoms, trajectory="ethanol.traj")
relax.run(fmax=0.001)

复制代码

运行结构优化
mpirun -n 16 gpaw python ethanol_relax.py共优化12步，29个SCF，用时4377s
运行结束后产生轨迹文件ethanol.traj和记录SCF过程的ethanol_relax.txt

计算能量差
完成优化后，读取优化完的结构计算气相与液相能量之差：

from ase.io import read as ase_read
from ase.io import write as ase_write
from gpaw import GPAW, PW
from ase.optimize import BFGSLineSearch
#from ase.filters import UnitCellFilter
from dftd3.ase import DFTD3
from ase.calculators.mixing import SumCalculator
from ase.constraints import FixAtoms
from ase.units import mol, kJ, kcal, Pascal, m
from ase.data.vdw import vdw_radii
from ase.build import molecule
from ase.parallel import parprint
from gpaw.solvation import (
# calculator
SolvationGPAW,
# cavities
EffectivePotentialCavity,
ADM12SmoothStepCavity,
#FG02SmoothStepCavity,
# custom classes for the cavities
Power12Potential,
ElDensity,
#SSS09Density,
# dielectric
LinearDielectric,
# non-el interactions
SurfaceInteraction,
VolumeInteraction,
LeakedDensityInteraction,
# surface and volume calculators
GradientSurface,
KB51Volume,
)
from ase.io.trajectory import Trajectory
atoms = Trajectory('ethanol.traj')[-1] #取最后一帧
# non-solvent related DFT parameters
h = 0.1
# solvent parameters for water from J. Chem. Phys. 141, 174108 (2014)
u0 = 0.180 # eV
epsinf = 78.36 # dimensionless
gamma = 18.4 * 1e-3 * Pascal * m # convert from dyne / cm to eV / Angstrom**2
T = 298.15 # Kelvin
p = -0.35 * 1e9 * Pascal
vdw_radii = vdw_radii.copy()
vdw_radii[1] = 1.09
vdw_radii[25] = 2.04
kappa_T = 4.53e-10 / Pascal
def atomic_radii(atoms_vdwr):
return [vdw_radii[n] for n in atoms_vdwr.numbers]
#液相
dft1 = SolvationGPAW(mode='lcao', basis='dzp', xc='PBE',nbands='110%',kpts=(2,2,2),
parallel=dict(band=2, # band parallelization
augment_grids=True, # use all cores for XC/Poisson
sl_auto=True, # enable ScaLAPACK parallelization
use_elpa=True),
h=h, txt='ethanol_water.txt',
# occupations={'name': 'fermi-dirac', 'width': 0.05},
cavity=EffectivePotentialCavity(
effective_potential=Power12Potential(atomic_radii, u0),
temperature=T,
surface_calculator=GradientSurface(),
volume_calculator=KB51Volume(compressibility=kappa_T,temperature=T)),
dielectric=LinearDielectric(epsinf=epsinf),
interactions=[SurfaceInteraction(surface_tension=gamma),
VolumeInteraction(pressure=p),
LeakedDensityInteraction(voltage=1)])
atoms.calc = SumCalculator([DFTD3(method="PBE", damping="d3bj"), dft1])
Ewater = atoms.get_potential_energy()
#气相
dft2=GPAW(xc='PBE',mode='lcao',basis='dzp',h=h,nbands='110%',kpts=(2,2,2),
parallel=dict(band=2, # band parallelization
augment_grids=True, # use all cores for XC/Poisson
sl_auto=True, # enable ScaLAPACK parallelization
use_elpa=True),
txt='ethanol_gas.txt')
atoms.calc = SumCalculator([DFTD3(method="PBE", damping="d3bj"), dft2])
Egasphase = atoms.get_potential_energy()
DGSol_eV = Ewater - Egasphase
DGSol_kJ_per_mol = DGSol_eV / (kJ / mol)
DGSol_kcal_per_mol = DGSol_eV / (kcal / mol)
parprint('calculated Delta Gsol = %.0f meV = %.1f kJ / mol = %.1f kcal / mol' %
(DGSol_eV * 1000., DGSol_kJ_per_mol, DGSol_kcal_per_mol))

复制代码

运行上述.py文件得到的结果为-0.185eV,与其他软件的结果相近。气相每次scf迭代用时1s，液相每次迭代用时10s。

ABACUS
安装
与GPAW一样，ABACUS也可以直接用conda安装所有需要的库与本体，建议使用新环境隔离开。
安装：conda create -n abacus_env abacus "libblas=*=*mkl" mpich -c conda-forge
运行：mpirun -n 16 abacus

模型构建与输入文件
ABACUS的结构文件STRU文件有点恶心，使用ABACUS的ASE插件来生成最方便，参考：https://gitlab.com/1041176461/ase-abacus。直接使用上一部分产生的初始.cif文件即可：

from ase.io import read as ase_read
from ase.io import write as ase_write
from ase.constraints import FixAtoms
from ase.build import molecule
ethanol = ase_read('ethanol_initial.cif')
#实际上ASE atom对象的FixAtoms和定义磁矩也是有效的，比如固定slab中的底层可以：
#mno2.set_constraint(FixAtoms(mask=mno2.positions[:, 2] < 3))
#表明固定二氧化锰z坐标低于3埃的原子
#磁矩，给Mn原子一个3的初猜：
#maganetic = []
#for atom in mno2:
# if atom.symbol == "Mn":
# maganetic.append(3)
# else:
# maganetic.append(0)
#mno2.set_initial_magnetic_moments(maganetic)
#ABACUS使用LCAO时需要同时提供赝势文件与数值轨道基组，可以在官网下载
pp={'C':'C_ONCV_PBE-1.0.upf',
'H':'H_ONCV_PBE-1.0.upf',
'O':'O_ONCV_PBE-1.0.upf'}
basis={'C':'C_gga_7au_100Ry_2s2p1d.orb',
'H':'H_gga_6au_100Ry_2s1p.orb',
'O':'O_gga_7au_100Ry_2s2p1d.orb'}
ase_write('STRU',ethanol,format='abacus',pp=pp,basis=basis)

复制代码

然后得到相应的STRU文件：

ATOMIC_SPECIES
C 12.011 C_ONCV_PBE-1.0.upf
O 15.999 O_ONCV_PBE-1.0.upf
H 1.008 H_ONCV_PBE-1.0.upf
NUMERICAL_ORBITAL
C_gga_7au_100Ry_2s2p1d.orb
O_gga_7au_100Ry_2s2p1d.orb
H_gga_6au_100Ry_2s1p.orb
LATTICE_CONSTANT
1.8897261258369282
LATTICE_VECTORS
14.062514000 0.0000000000 0.0000000000
0.0000000000 12.244742000 0.0000000000
0.0000000000 0.0000000000 11.773866000
ATOMIC_POSITIONS
Direct
C
0.0000000000
2
0.5770520000 0.4603490000 0.5000000000 1 1 1 mag 0.0
0.4939820000 0.5387370000 0.5000000000 1 1 1 mag 0.0
O
0.0000000000
1
0.4093540000 0.4744540000 0.5000000000 1 1 1 mag 0.0
H
0.0000000000
6
0.3555550000 0.5242050000 0.5000000000 1 1 1 mag 0.0
0.4970080000 0.5916610000 0.5753310000 1 1 1 mag 0.0
0.4970080000 0.5916610000 0.4246690000 1 1 1 mag 0.0
0.6444450000 0.5048730000 0.5000000000 1 1 1 mag 0.0
0.5742370000 0.4083390000 0.5752410000 1 1 1 mag 0.0
0.5742370000 0.4083390000 0.4247590000 1 1 1 mag 0.0

复制代码

INPUT文件:

INPUT_PARAMETERS
suffix ethanol
ntype 3
nelec 0.0
nspin 1
pseudo_dir /path_to_pp/ #赝势文件夹
orbital_dir /path_to_orb/ #轨道基组文件夹
dft_functional pbe
ecutwfc 100 # Rydberg #LCAO使用建议值100Ry,使用PW方法最好测试下
scf_thr 1e-7 # Rydberg #用比较高的收敛限
basis_type lcao #使用LCAO
calculation relax # this is the key parameter telling abacus to do a scf calculation
force_thr_ev 0.01 # the threshold of the force convergence, in unit of eV/Angstrom
relax_nmax 100 # the maximal number of ionic iteration steps
out_stru 1 #输出结构
vdw_method d3_bj #dftd3-bj校正
#溶剂模型参数
imp_sol 1 #打开溶剂模型
eb_k 80
tau 0.000010798
sigma_k 0.6
nc_k 0.00037

复制代码

KPT文件:

K_POINTS
0
Gamma
2 2 2 0 0 0

复制代码

和VASP类似（幸好不用合并赝势文件），把他们放到同一个文件夹，运行mpirun -n 16 abacus，得到最后一帧结构。

共优化27步，用时7785s

溶解自由能
使用优化后的结构，分别计算气相与液相的能量。

INPUT_PARAMETERS
suffix ethanol
ntype 3
nelec 0.0
nspin 1
pseudo_dir /path_to_pp/ #赝势文件夹
orbital_dir /path_to_orb/ #轨道基组文件夹
dft_functional pbe
ecutwfc 100 # Rydberg #LCAO使用建议值100Ry,使用PW方法最好测试下
scf_thr 1e-7 # Rydberg #用比较高的收敛限
basis_type lcao #使用LCAO
calculation scf # this is the key parameter telling abacus to do a scf calculation
vdw_method d3_bj #dftd3-bj校正
#溶剂模型参数
imp_sol 0 #关闭溶剂模型
eb_k 80
tau 0.000010798
sigma_k 0.6
nc_k 0.00037

复制代码

气相中每次scf迭代用时约为0.8s，总能量-841.1220375167616 eV 。液相中每次scf迭代用时16s，总能量-841.3260286872329 eV，溶解自由能：-0.20 eV符合。

QE
安装
QE要使用溶剂模型需要在编译时加入Environ插件。由于是AMD平台，不考虑的MKL的话编译很简单，先进入解压后的Environ源代码文件夹，输入：./configure --with-qe=<absolute_path_to_QE>，成功后进行编译：make -j compile。
没有报错后进入QE的目录，输入./configure --with-environ=<absolute_path_to_Environ> -enable-openmp，成功后开始编译QEmake -jN all。详细QE编译并且加上MKL参考http://sobereva.com/562。在测试中发现有些版本QE的configure无法正确识别MKL，如果出现此情况，有条件的可以自行修改，没条件的又不急需那个版本的话直接换一个版本的源码包。

结构优化(qe用的人比较多，就不放输入文件了)
使用SSSP Efficiency的赝势库，截断能设置为60，600。
Environ input （environ.in）:

&ENVIRON
verbose = 1
environ_type = "water" #应用默认的参数
environ_thr = 1.d2
/
&BOUNDARY
solvent_mode = 'full'
/
&ELECTROSTATIC
pbc_dim = 3

复制代码

将.in文件和environ.in放到一个文件夹运行：

mpirun -n 16 pw.x < enthanol.in -environ液相优化用时494s，12步。每次scf迭代8s。气相每次迭代3.8秒。

溶解自由能
液相中总能量-85.6621819080 Ry ，气相中-85.65146975 Ry，溶解自由能为-0.14 eV，差距有点多。

CP2K
安装
参考http://sobereva.com/586
直接使用Multiwfn生成的输入文件默认值进行优化。首先使用非周期性，gamma点进行计算。中间偶尔几步能量波动特别大，在优化85步之后，SCF直接由于solvation energy异常（+1022 kcal/mol）不收敛了。
现在改成三维周期性，直接用Multiwfn生成的输入文件（截断能：55，400；Perodic求解器；RHO_MIN 0.0001841； RHO_MAX 0.0013604）。同样，优化前几步就出现很大的能量波动以及不合理的solvation energy：

将输入文件中的值改成和Envrion中的值一致(Dupont et al., 2013)，即α+γ 50，β -0.35，RHO_MIN 0.0001，RHO_MAX 0.005。将截断能提高至60、800再次进行结构优化。这次经过8步后优化收敛，每步scf迭代为10s左右，总用时1744s，Solvation energy为-0.17 eV符合。

小结
以上软件应该都能够通过调整参数来得到合理的小分子自由能，但是cp2k需要很小心地测试隐式溶剂模型截断能（测试力？），有时候虽然第一步能给个负值，但是在优化过程中就会不收敛，尽管此时构象并没有发生明显变化。其他软件在这一环节中没有出现cp2k的问题。优化后GPAW给出的结构与pyscf优化的差距稍大，部分键长差距2%左右，其他软件优化的差距在1%以内。

三、slab吸附模型测试
参考文献中优化好的晶胞(Boyd et al., 2021)，进行扩胞搞了个很大的双层MnO2 slab，原子数448，层间有一些水和钾离子。
1. GPAW
GPAW目前的隐式溶剂模型不支持带电体系；偶极校正需要校正方向与slab面垂直且在该方向上为非周期性边界；带电体系不支持周期性和非周期性边界的混合。因此，目前的GPAW模拟水溶液中的吸附、带电体系的吸附有很大的局限。
2. ABACUS
正在施工。。。。有空补上
3. QE
就不难为我小小的机子了，直接放弃
4. CP2K
用OT很勉强，截断能不能很高。。。。。有空补上

		自动登录 Automatic login	找回密码 Forget password
密码 Password			注册 Register

[综合讨论] 免费第一性原理软件的溶剂模型测试笔记

评分 Rate