Dear ASE users and developers:
I am new to run CP2K coupled with ASE, here I learn ASE communicate with cp2k_shell via stdin/stdout. The example scripts I found on the internet set the other sections in the cp2k calculator and leave the SUBSYS section empty for the ASE atoms class.
Here I want to know is it possible to power DFT+U calculations in CP2K by ASE? Because the configurations of the DFT+U method and the important BS section are both within the SUBSYS/KIND section, and the KIND secation also allows one to define multiple atom kinds for one kind of element, the ASE atoms instance don’t have these properties.
Any reply will be helpful!
Yours,
Youmu
Dear all,
After diving into the source code of ase.calculators.cp2k, I found the only limitation here is the way CP2K class sets KIND section does not match CP2K program well.
In CP2K, the keyword of KIND section can be any string and it is the subkeyword ELEMENT that determines the real element type. For example, in \mathrm{Fe_3O_4}, we get three different kinds of \mathrm{Fe} and each of them has different atomic configuration, which can be set explicitly in KIND/BS section. However, the ase.Atom class does not allow to arbitarily set the ase.Atom.symbol because ase gets chemical symbols through a hash map between symbols and atomic numbers.
So my solution is to add a new property kind (kinds) to the ase.Atom (ase.Atoms), respectively. The kind of one atom instance can be set as any string and its default value its chemical symbol. Furthermore, an Atoms instance can set (or get) its kinds through Atoms.set_kinds (or Atoms.get_kinds) method. Having these properties added to ase, the only obstacle is in ase.calculators.CP2K._generate_input which is about writing KIND section and coords in the generated .inp file. I replaced the self.atoms.get_chemical_symbols() to the new get_kinds() method and added a new keyword ELEMENT. Also the elements in writting coords are replaced by the kinds.
Thanks to the well formatted ase codes, I was able to add features I need. After doing these “hacking” stuffs, I provide a working example below, containing a Python script and the generated CP2K inp file.
Since this modification has not been well tested, it will take some time before I open it (or being asked). If anyone is interested in this feature, please contact me.
Best regards,
Youmu
Script:
#!~/anaconda3/envs/ASE/bin/python
from ase.calculators.cp2k import CP2K
from ase.build import molecule
calc = CP2K(directory='H2O')
atoms = molecule('H2O')
atoms.set_kinds(['O', 'H1', 'H2'])
atoms.calc = calc
print(atoms.get_kinds())
atoms.center(vacuum=2.0)
print(atoms.get_potential_energy())
Generated input:
&GLOBAL
PROJECT H2O/cp2k
PRINT_LEVEL LOW
&END GLOBAL
&FORCE_EVAL
METHOD Quickstep
STRESS_TENSOR ANALYTICAL
&PRINT
&STRESS_TENSOR ON
&END STRESS_TENSOR
&END PRINT
&DFT
BASIS_SET_FILE_NAME BASIS_MOLOPT
POTENTIAL_FILE_NAME POTENTIAL
&MGRID
CUTOFF [eV] 5.442277204873448682e+03
&END MGRID
&SCF
MAX_SCF 50
&END SCF
&LS_SCF
MAX_SCF 50
&END LS_SCF
&XC
&XC_FUNCTIONAL
&PADE
&END PADE
&END XC_FUNCTIONAL
&END XC
&POISSON
PERIODIC NONE
PSOLVER MT
&END POISSON
&END DFT
&SUBSYS
&COORD
O 2.000000000000000000e+00 2.763239000000000445e+00 2.596308999999999756e+00
H1 2.000000000000000000e+00 3.526478000000000446e+00 1.999999999999999778e+00
H2 2.000000000000000000e+00 2.000000000000000444e+00 1.999999999999999778e+00
&END COORD
&CELL
PERIODIC NONE
A 4.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00
B 0.000000000000000000e+00 5.526478000000000002e+00 0.000000000000000000e+00
C 0.000000000000000000e+00 0.000000000000000000e+00 4.596308999999999756e+00
&END CELL
&KIND O
ELEMENT O
BASIS_SET DZVP-MOLOPT-SR-GTH
POTENTIAL GTH-LDA
&END KIND
&KIND H1
ELEMENT H
BASIS_SET DZVP-MOLOPT-SR-GTH
POTENTIAL GTH-LDA
&END KIND
&KIND H2
ELEMENT H
BASIS_SET DZVP-MOLOPT-SR-GTH
POTENTIAL GTH-LDA
&END KIND
&END SUBSYS
&END FORCE_EVAL