Unexpected error when using Crystallographic Information Framework (CIF) File Input

Question

Unexpected error when using Crystallographic Information Framework (CIF) File Input

Closed this issue a year ago · 5 comments

Disclaimer

Importing a Quantum Espresso Input works perfectly fine

Prerequisite

I have a CIF file of Graphene (P6/mmm) [unit cell image like above] in the same folder as the following which I want to import in my DFTK code to run-

data_Graphene
_audit_creation_date              2023-06-30
_audit_creation_method            'Materials Studio'
_symmetry_space_group_name_H-M    'P1'
_symmetry_Int_Tables_number       1
_symmetry_cell_setting            triclinic
loop_
_symmetry_equiv_pos_as_xyz
  x,y,z
_cell_length_a                    2.4595
_cell_length_b                    4.2600
_cell_length_c                    30.0000
_cell_angle_alpha                 90.0000
_cell_angle_beta                  90.0000
_cell_angle_gamma                 90.0000
loop_
_atom_site_label
_atom_site_type_symbol
_atom_site_fract_x
_atom_site_fract_y
_atom_site_fract_z
_atom_site_U_iso_or_equiv
_atom_site_adp_type
_atom_site_occupancy
C1     C     0.00000   0.33333  -0.00000   0.00000  Uiso   1.00
C2     C     0.50000   0.83333   0.00000   0.00000  Uiso   1.00
C3     C    -0.00000   0.66667   0.00000   0.00000  Uiso   1.00
C4     C     0.50000   0.16667   0.00000   0.00000  Uiso   1.00
loop_
_geom_bond_atom_site_label_1
_geom_bond_atom_site_label_2
_geom_bond_distance
_geom_bond_site_symmetry_2
_ccdc_geom_bond_type
C1     C3      1.420   .     A
C1     C4      1.420   1_455 A
C2     C4      1.420   1_565 A
C2     C3      1.420   .     A
C2     C3      1.420   1_655 A
C3     C2      1.420   1_455 A
C4     C2      1.420   1_545 A
C4     C1      1.420   1_655 A

Problem:

First Try

I was trying to import a CIF file of Graphene with the help of documentation. So, I wrote the following code-

using DFTK, AtomsIO, Plots, Unitful, UnitfulAtomic
system = load_system("Graphene.cif")
system = attach_psp(system, C="hgh/pbe/c-q4.hgh", )
model = model_PBE(system; temperature=0.01)
basis = PlaneWaveBasis(model; Ecut=30, kgrid=(2, 2, 2))
ρ0 = guess_density(basis, system)
scfres = self_consistent_field(basis, ρ=ρ0);

The above code gives the following error for line number 3-

No pseudo identifier given for element C1.

Stacktrace:
 [1] error(s::String)
   @ Base [.\error.jl:35](https://file+.vscode-resource.vscode-cdn.net/d%3A/Programming/Julia%20Workshop/My%20Code/error.jl:35)
 [2] (::DFTK.var"#754#755"{Dict{Symbol, String}})(atom::Atom{3, Quantity{Float64, 𝐋, Unitful.FreeUnits{(Å,), 𝐋, nothing}}, Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(Å, ps^-1), 𝐋 𝐓^-1, nothing}}, Quantity{Float64, 𝐌, Unitful.FreeUnits{(u,), 𝐌, nothing}}})
   @ DFTK [C:\Users\Aritra\.julia\packages\DFTK\k0wtY\src\pseudo\attach_psp.jl:33](file:///C:/Users/Aritra/.julia/packages/DFTK/k0wtY/src/pseudo/attach_psp.jl:33)
 [3] iterate
   @ [.\generator.jl:47](https://file+.vscode-resource.vscode-cdn.net/d%3A/Programming/Julia%20Workshop/My%20Code/generator.jl:47) [inlined]
 [4] collect(itr::Base.Generator{FlexibleSystem{3, Atom{3, Quantity{Float64, 𝐋, Unitful.FreeUnits{(Å,), 𝐋, nothing}}, Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(Å, ps^-1), 𝐋 𝐓^-1, nothing}}, Quantity{Float64, 𝐌, Unitful.FreeUnits{(u,), 𝐌, nothing}}}, Quantity{Float64, 𝐋, Unitful.FreeUnits{(Å,), 𝐋, nothing}}}, DFTK.var"#754#755"{Dict{Symbol, String}}})
   @ Base [.\array.jl:782](https://file+.vscode-resource.vscode-cdn.net/d%3A/Programming/Julia%20Workshop/My%20Code/array.jl:782)
 [5] map
   @ [.\abstractarray.jl:3289](https://file+.vscode-resource.vscode-cdn.net/d%3A/Programming/Julia%20Workshop/My%20Code/abstractarray.jl:3289) [inlined]
 [6] attach_psp(system::FlexibleSystem{3, Atom{3, Quantity{Float64, 𝐋, Unitful.FreeUnits{(Å,), 𝐋, nothing}}, Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(Å, ps^-1), 𝐋 𝐓^-1, nothing}}, Quantity{Float64, 𝐌, Unitful.FreeUnits{(u,), 𝐌, nothing}}}, Quantity{Float64, 𝐋, Unitful.FreeUnits{(Å,), 𝐋, nothing}}}, pspmap::Dict{Symbol, String})
   @ DFTK [C:\Users\Aritra\.julia\packages\DFTK\k0wtY\src\pseudo\attach_psp.jl:26](file:///C:/Users/Aritra/.julia/packages/DFTK/k0wtY/src/pseudo/attach_psp.jl:26)
 [7] #attach_psp#756
   @ [C:\Users\Aritra\.julia\packages\DFTK\k0wtY\src\pseudo\attach_psp.jl:41](file:///C:/Users/Aritra/.julia/packages/DFTK/k0wtY/src/pseudo/attach_psp.jl:41) [inlined]
 [8] top-level scope
   @ [d:\Programming\Julia](file:///D:/Programming/Julia) Workshop\My Code\dft-graphene.ipynb:1

Second Try

As the previous code said, No pseudo identifier given for element C1, I made some changed in line 3 like following-

using DFTK, AtomsIO, Plots, Unitful, UnitfulAtomic
system = load_system("Graphene.cif")
system = attach_psp(system, C1="hgh/pbe/c-q4.hgh", C2="hgh/pbe/c-q4.hgh", C3="hgh/pbe/c-q4.hgh", C4="hgh/pbe/c-q4.hgh")
model = model_PBE(system; temperature=0.01)
basis = PlaneWaveBasis(model; Ecut=30, kgrid=(2, 2, 2))
ρ0 = guess_density(basis, system)
scfres = self_consistent_field(basis, ρ=ρ0);

This updated version passes the 3rd line, but throws an error at line 4-

KeyError: key :C1 not found

Stacktrace:
  [1] getindex(h::Dict{Symbol, PeriodicTable.Element}, key::Symbol)
    @ Base [.\dict.jl:484](https://file+.vscode-resource.vscode-cdn.net/d%3A/Programming/Julia%20Workshop/My%20Code/dict.jl:484)
  [2] getindex(e::PeriodicTable.Elements, i::Symbol)
    @ PeriodicTable [C:\Users\Aritra\.julia\packages\PeriodicTable\aQMqe\src\PeriodicTable.jl:160](file:///C:/Users/Aritra/.julia/packages/PeriodicTable/aQMqe/src/PeriodicTable.jl:160)
  [3] ElementPsp(key::Symbol; psp::PspHgh{Float64})
    @ DFTK [C:\Users\Aritra\.julia\packages\DFTK\k0wtY\src\elements.jl:94](file:///C:/Users/Aritra/.julia/packages/DFTK/k0wtY/src/elements.jl:94)
  [4] (::DFTK.var"#758#764"{Atom{3, Quantity{Float64, 𝐋, Unitful.FreeUnits{(Å,), 𝐋, nothing}}, Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(Å, ps^-1), 𝐋 𝐓^-1, nothing}}, Quantity{Float64, 𝐌, Unitful.FreeUnits{(u,), 𝐌, nothing}}}, String})()
    @ DFTK [C:\Users\Aritra\.julia\packages\DFTK\k0wtY\src\external\atomsbase.jl:25](file:///C:/Users/Aritra/.julia/packages/DFTK/k0wtY/src/external/atomsbase.jl:25)
  [5] get!(default::DFTK.var"#758#764"{Atom{3, Quantity{Float64, 𝐋, Unitful.FreeUnits{(Å,), 𝐋, nothing}}, Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(Å, ps^-1), 𝐋 𝐓^-1, nothing}}, Quantity{Float64, 𝐌, Unitful.FreeUnits{(u,), 𝐌, nothing}}}, String}, h::Dict{String, ElementPsp}, key::String)
    @ Base [.\dict.jl:468](https://file+.vscode-resource.vscode-cdn.net/d%3A/Programming/Julia%20Workshop/My%20Code/dict.jl:468)
  [6] (::DFTK.var"#757#763"{Dict{String, ElementPsp}})(atom::Atom{3, Quantity{Float64, 𝐋, Unitful.FreeUnits{(Å,), 𝐋, nothing}}, Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(Å, ps^-1), 𝐋 𝐓^-1, nothing}}, Quantity{Float64, 𝐌, Unitful.FreeUnits{(u,), 𝐌, nothing}}})
    @ DFTK [C:\Users\Aritra\.julia\packages\DFTK\k0wtY\src\external\atomsbase.jl:24](file:///C:/Users/Aritra/.julia/packages/DFTK/k0wtY/src/external/atomsbase.jl:24)
  [7] iterate
    @ [.\generator.jl:47](https://file+.vscode-resource.vscode-cdn.net/d%3A/Programming/Julia%20Workshop/My%20Code/generator.jl:47) [inlined]
  [8] collect(itr::Base.Generator{FlexibleSystem{3, Atom{3, Quantity{Float64, 𝐋, Unitful.FreeUnits{(Å,), 𝐋, nothing}}, Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(Å, ps^-1), 𝐋 𝐓^-1, nothing}}, Quantity{Float64, 𝐌, Unitful.FreeUnits{(u,), 𝐌, nothing}}}, Quantity{Float64, 𝐋, Unitful.FreeUnits{(Å,), 𝐋, nothing}}}, DFTK.var"#757#763"{Dict{String, ElementPsp}}})
    @ Base [.\array.jl:782](https://file+.vscode-resource.vscode-cdn.net/d%3A/Programming/Julia%20Workshop/My%20Code/array.jl:782)
  [9] map
    @ [.\abstractarray.jl:3289](https://file+.vscode-resource.vscode-cdn.net/d%3A/Programming/Julia%20Workshop/My%20Code/abstractarray.jl:3289) [inlined]
 [10] parse_system(system::FlexibleSystem{3, Atom{3, Quantity{Float64, 𝐋, Unitful.FreeUnits{(Å,), 𝐋, nothing}}, Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(Å, ps^-1), 𝐋 𝐓^-1, nothing}}, Quantity{Float64, 𝐌, Unitful.FreeUnits{(u,), 𝐌, nothing}}}, Quantity{Float64, 𝐋, Unitful.FreeUnits{(Å,), 𝐋, nothing}}})
    @ DFTK [C:\Users\Aritra\.julia\packages\DFTK\k0wtY\src\external\atomsbase.jl:19](file:///C:/Users/Aritra/.julia/packages/DFTK/k0wtY/src/external/atomsbase.jl:19)
 [11] model_PBE(::FlexibleSystem{3, Atom{3, Quantity{Float64, 𝐋, Unitful.FreeUnits{(Å,), 𝐋, nothing}}, Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(Å, ps^-1), 𝐋 𝐓^-1, nothing}}, Quantity{Float64, 𝐌, Unitful.FreeUnits{(u,), 𝐌, nothing}}}, Quantity{Float64, 𝐋, Unitful.FreeUnits{(Å,), 𝐋, nothing}}}; kwargs::Base.Pairs{Symbol, Float64, Tuple{Symbol}, NamedTuple{(:temperature,), Tuple{Float64}}})
    @ DFTK [C:\Users\Aritra\.julia\packages\DFTK\k0wtY\src\standard_models.jl:81](file:///C:/Users/Aritra/.julia/packages/DFTK/k0wtY/src/standard_models.jl:81)
 [12] top-level scope
    @ [d:\Programming\Julia](file:///D:/Programming/Julia) Workshop\My Code\dft-graphene.ipynb:1

Third Try

Following the recent error I made some changed in line 3 like following (there is nothing to update in line 4)-

using DFTK, AtomsIO, Plots, Unitful, UnitfulAtomic
system = load_system("Graphene.cif")
system = attach_psp(system, :C1="hgh/pbe/c-q4.hgh", :C2="hgh/pbe/c-q4.hgh", :C3="hgh/pbe/c-q4.hgh", :C4="hgh/pbe/c-q4.hgh")
model = model_PBE(system; temperature=0.01)
basis = PlaneWaveBasis(model; Ecut=30, kgrid=(2, 2, 2))
ρ0 = guess_density(basis, system)
scfres = self_consistent_field(basis, ρ=ρ0);

In this case, I got a syntax error like below-

syntax: invalid keyword argument name ":C1" around [d:\Programming\Julia](file:///D:/Programming/Julia) Workshop\My Code\dft-graphene.ipynb:1

Stacktrace:
 [1] top-level scope
   @ [d:\Programming\Julia](file:///D:/Programming/Julia) Workshop\My Code\dft-graphene.ipynb:1

Changing Approach & Solution:

I tried changing the CIF file such as from C1, C2, C3 and C4 to C which looks like following-

data_Graphene
_audit_creation_date              2023-06-30
_audit_creation_method            'Materials Studio'
_symmetry_space_group_name_H-M    'P1'
_symmetry_Int_Tables_number       1
_symmetry_cell_setting            triclinic
loop_
_symmetry_equiv_pos_as_xyz
  x,y,z
_cell_length_a                    2.4595
_cell_length_b                    4.2600
_cell_length_c                    30.0000
_cell_angle_alpha                 90.0000
_cell_angle_beta                  90.0000
_cell_angle_gamma                 90.0000
loop_
_atom_site_label
_atom_site_type_symbol
_atom_site_fract_x
_atom_site_fract_y
_atom_site_fract_z
_atom_site_U_iso_or_equiv
_atom_site_adp_type
_atom_site_occupancy
C     C     0.00000   0.33333  -0.00000   0.00000  Uiso   1.00
C     C     0.50000   0.83333   0.00000   0.00000  Uiso   1.00
C     C    -0.00000   0.66667   0.00000   0.00000  Uiso   1.00
C     C     0.50000   0.16667   0.00000   0.00000  Uiso   1.00
loop_
_geom_bond_atom_site_label_1
_geom_bond_atom_site_label_2
_geom_bond_distance
_geom_bond_site_symmetry_2
_ccdc_geom_bond_type
C     C      1.420   .     A
C     C      1.420   1_455 A
C     C      1.420   1_565 A
C     C      1.420   .     A
C     C      1.420   1_655 A
C     C      1.420   1_455 A
C     C      1.420   1_545 A
C     C      1.420   1_655 A

This worked smoothly and gave the expected results.

Trying Other Structures for Confirmation

When I found the possible solution, I tried with some other structures from Crystallography Open Database
to confirm the solution.

No.	Elements	CIF Download Link	Used Code (3rd Line: Attaching Pseudopotentials)	Encountered Error	Status After Implementing This Solution (Editing CIF File)
1	Ba O4 S	🔗	system = attach_psp(system, Ba="hgh/pbe/ba-q10.hgh", S="hgh/pbe/s-q6.hgh", O="hgh/pbe/o-q6.hgh")	No pseudo identifier given for element Ba1.	✅
2	H5 O7 P V	🔗	system = attach_psp(system, H="hgh/pbe/h-q1.hgh", O="hgh/pbe/o-q6.hgh", P="hgh/pbe/p-q5.hgh", V="hgh/pbe/v-q13.hgh")	No pseudo identifier given for element V1.	✅

Maybe a PR?

I hope this will help both computational and experimental researchers as CIF is the most common file used by both of them when it comes to materials science. So, I can make a PR by adding a suitable example to the documentation regarding this.

Answer 1 · 2023-07-02T09:42:13.000Z

Thanks @aritraroy24 for this detailed report. User experiences like this are highly appreciated to make it easier for everyone to work with DFTK.

In my opinion what you try to do here should just work without manual editing. The reason it does not boils down to the following:
AtomsIO.jl (which provides load_system) uses Chemfiles.jl to parse the Cif file. This passes on the bare labels c1 etc without stripping off the tailling numbers, such that dftk thinks the element symbol is c1.

There are multiple ways to fix this:

in chemfiles
in AtomsIO
in DFTK

To give us another datapoint, could you check parsing works with the original files and the ASE parser (load_system(ASE(), filename)) - check AtomsIO docs for details. You will need the AtomsIOPython package for this.

Answer 2 · 2023-07-02T09:48:43.000Z

Could you also check if #864 fixes the issue?

Answer 3 · 2023-07-02T11:14:38.000Z

I think this is actually an issue on the DFTK side. We use the atomic_symbol, where one is supposed to use the atomic_number to identify an element. This is fixed in #864. I have now also added a test and I will merge it once this runs through.

Answer 4 · 2023-07-02T15:01:30.000Z

@mfherbst, Happy to report the bug and thanks for the quick action. I'm actually out of the station and will be available from Tuesday. I'll try #864 on Tuesday and will give the update.

Answer 5 · 2023-07-02T15:23:09.000Z

Sounds good. Feel free to reopen this in case my PR has not fixed the bug.