Merge pull request #5 from ejmeitz/ucposcar

add IO for uc and ssposcar

Author: ejmeitz

Project.toml

@@ -1,19 +1,21 @@
 name = "SimpleCrystals"
 uuid = "64031d72-e220-11ed-1a7e-43a2532b2fa8"
 authors = ["Ethan Meitz <emeitz@andrew.cmu.edu>"]
-version = "0.4.0"
+version = "0.4.1"
 
 [deps]
 AtomsBase = "a963bdd2-2df7-4f54-a1ee-49d51e6be12a"
 PeriodicTable = "7b2266bf-644c-5ea3-82d8-af4bbd25a884"
+Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
 [compat]
 AtomsBase = "0.5"
 PeriodicTable = "1.1"
+Printf = "1.11.0"
 Reexport = "1"
 StaticArrays = "1.5.17"
 Unitful = "1"
-julia = "1.7,1.10"
+julia = "1.10"

README.md

@@ -140,3 +140,6 @@ for atom in fcc_crystal
     #do something
 end
 ```
+
+There are also functions `to_ucposcar` and `to_ssposcar` to generate poscar files. If your crystal
+does not use Unitful ensure your positions and lattice vectors are in Angstroms. 

src/FileWriter.jl

@@ -1,4 +1,4 @@
-export to_xyz
+export to_xyz, to_ucposcar, to_ssposcar
 
 function to_xyz(crystal::Crystal{3}, outpath)
     N_atoms = length(crystal)
@@ -22,4 +22,94 @@ function to_xyz(crystal::Crystal{2}, outpath)
         println(file,"$(string(atomic_symbol(atom))) $(ustrip(atom.position[1])) $(ustrip(atom.position[2]))")
     end
     close(file)
+end
+
+"""
+    cartesian_to_fractional(latt::BravaisLattice{3}, x::AbstractVector) where T
+
+Convert a Cartesian coordinate `x` into fractional coordinates w.r.t. the
+primitive vectors of `latt`.
+"""
+cartesian_to_fractional(lattice_vectors::AbstractMatrix, x::AbstractVector) = lattice_vectors \ x
+
+
+function convert_to_ang_and_strip(x)
+    if eltype(x) <: Unitful.Length
+        x = ustrip.(u"Å", x)
+    else
+        @error "Expected length units, got $(unit(first(x)))"
+    end
+end
+
+function write_poscar(cryst::Crystal{D}, cell::AbstractMatrix, path::AbstractString) where D
+    # 1) collect element order & counts
+    all_syms = [atomic_symbol(atom) for atom in cryst.atoms]
+    # preserve order of first appearance in basis
+    order = unique([atomic_symbol(atom) for atom in cryst.basis])
+    counts = [count(==(sym), all_syms) for sym in order]
+
+    open(path, "w") do io
+        println(io, "POSCAR")                   # comment line
+        println(io, "1.0")                      # universal scale
+
+        # primitive vectors as rows
+        for i in 1:D
+            v = convert_to_ang_and_strip(cell[:,i])
+            for el in v
+                @printf(io, "%.15f ", el)
+            end
+            println(io)
+        end
+
+        # elements and counts
+        println(io, join(order, " "))
+        println(io, join(counts, " "))
+
+        println(io, "Direct")  # fractional coords
+
+        # write atoms in block order
+        for sym in order
+            for atom in cryst.atoms
+                if atomic_symbol(atom) == sym
+                    x = convert_to_ang_and_strip(position(atom))
+                    vec = convert_to_ang_and_strip(cell)
+                    x_frac = cartesian_to_fractional(vec, x)
+                    for el in x_frac
+                        @printf(io, "%.15f ", el)
+                    end
+                    println(io)
+                end
+            end
+        end
+    end
+end
+
+"""
+to_ucposcar(cryst::Crystal{3}; filename="ucposcar")
+
+Automatically extracts the *unit cell* from `cryst` (i.e. N_unit_cells ≡ (1,1,1))
+and writes it to a POSCAR-style file named `filename`.
+"""
+function to_ucposcar(cryst::Crystal{D}, filename::AbstractString="ucposcar") where D
+    
+    n_cells = @SVector ones(Int, D)
+
+    uc = Crystal(
+        cryst.lattice,
+        cryst.basis,
+        n_cells,
+        cryst.basis
+    )
+    write_poscar(uc, cryst.lattice.primitive_vectors, filename)
+end
+
+"""
+    supercell(cryst::Crystal{3}; filename="ssposcar")
+
+Writes the full supercell defined by `cryst.N_unit_cells` to a
+POSCAR-style file named `filename`.
+"""
+function to_ssposcar(cryst::Crystal, filename::AbstractString="ssposcar")
+    cell = cryst.lattice.primitive_vectors .* cryst.N_unit_cells
+    write_poscar(cryst, cell, filename)
 end

src/SimpleCrystals.jl

@@ -2,6 +2,7 @@ module SimpleCrystals
 
 using Reexport
 using PeriodicTable
+using Printf
 
 @reexport using AtomsBase
 @reexport using Unitful