Skip to content

Commit

Permalink
Some modification and new files for QHA (#294)
Browse files Browse the repository at this point in the history 
* Add file for generating deformations for QHA.

* Update qha.py: Added tref and a new method for calculating the thermal expansion coefficient.

* Add a script for v-ZSISA-QHA . "Physical Review B 110 (1), 014103"

* modification in qha_aproximation.py.

* Some modification on deformation_utils and qha.

* Some modifications
  • Loading branch information
@samarerstm
samarerstm authored Oct 15, 2024
1 parent 2747d22 commit 567fa3f
Show file tree
Hide file tree
Showing 3 changed files with 2,390 additions and 6 deletions.
193 changes: 193 additions & 0 deletions abipy/dfpt/deformation_utils.py
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,193 @@
# deformation_utils.py

import numpy as np
from pymatgen.core import Structure, Lattice, Element
from abipy.core.symmetries import AbinitSpaceGroup
from pymatgen.symmetry.analyzer import SpacegroupAnalyzer
from abipy.abio.inputs import AbinitInput
import re


def generate_deformations_volumic(structure, eps_V=0.02, scales=None):
if scales is None:
scales = [-1, 0, 1, 2, 3]
rprim = structure.lattice.matrix
structures_new = {}

for i in scales:
rprim2 = np.copy(rprim)
rprim2[:, :] = rprim[:, :] * (1.00 + eps_V * i)**(1/3.)

structure2 = structure.copy()
structure2.lattice = Lattice(rprim2)
#structure2.scale_lattice(structure2.volume*(1.00 + eps_V * i))
namei = int(round(1000 * (1.00 + eps_V * i)))
formatted_namei = f"{namei:04d}"
structures_new[formatted_namei] = structure2


return structures_new

def generate_deformations(structure , eps=0.005):
spgrp = AbinitSpaceGroup.from_structure(structure )
print (spgrp)
spgrp_number=spgrp.spgid
rprim= structure.lattice.matrix

rprim2 = np.copy(rprim)
rprim_new = {}
structures_new = {}

if 1 <= spgrp_number <= 2:
disp=[[1,1,1,1,1,1], [0,1,1,1,1,1], [2,1,1,1,1,1], [1,0,1,1,1,1], [1,2,1,1,1,1], [1,1,0,1,1,1],
[1,1,2,1,1,1], [1,1,1,0,1,1], [1,1,1,2,1,1], [1,1,1,1,0,1], [1,1,1,1,2,1], [1,1,1,1,1,0],
[1,1,1,1,1,2], [0,0,1,1,1,1], [1,0,0,1,1,1], [1,1,0,0,1,1], [1,1,1,0,0,1], [1,1,1,1,0,0],
[0,1,0,1,1,1], [0,1,1,0,1,1], [0,1,1,1,0,1], [0,1,1,1,1,0], [1,0,1,0,1,1], [1,0,1,1,0,1],
[1,0,1,1,1,0], [1,1,0,1,0,1], [1,1,0,1,1,0], [1,1,1,0,1,0] , [0 ,0,0,0,0,0]]
if abs(rprim[1, 0]) > 1e-9 or abs(rprim[2, 0]) > 1e-9 or abs(rprim[2, 1]) > 1e-9:
print("Warning: The lattice is oriented such that xz =xy =yz =0 .")
rprim0 = np.copy(rprim)
a=rprim[0, :]
b=rprim[1, :]
c=rprim[2, :]
norm_a = np.linalg.norm(a)
norm_b = np.linalg.norm(b)
norm_c = np.linalg.norm(c)

# Compute angles between vectors
cos_ab = np.dot(a, b) / (norm_a * norm_b)
cos_ac = np.dot(a, c) / (norm_a * norm_c)
cos_bc = np.dot(b, c) / (norm_b * norm_c)

rprim0[0,0] = 1.0
rprim0[0,1] = 0.0
rprim0[0,2] = 0.0
rprim0[1,0] = cos_ab
rprim0[1,1] = np.sqrt(1-cos_ab**2)
rprim0[1,2] = 0.0
rprim0[2,0] = cos_ac
rprim0[2,1] = (cos_bc-rprim0[1,0]*rprim0[2,0])/rprim0[1,1]
rprim0[2,2] = np.sqrt(1.0-rprim0[2,0]**2-rprim0[2,1]**2)
rprim0[0,:] = rprim0[0,:]*norm_a
rprim0[1,:] = rprim0[1,:]*norm_b
rprim0[2,:] = rprim0[2,:]*norm_c
print("Old rprim:")
print(rprim)
print("New rprim:")
print(rprim0)

for pair in disp:
i,j,k,l,m,n = pair
rprim2[ :,0] = rprim0[ :,0] * (1.00 + eps * i) + rprim0[ :,1] * (eps * l) +rprim0[ :,2] * (eps * m)
rprim2[ :,1] = rprim0[ :,1] * (1.00 + eps * j) + rprim0[ :,2] * (eps * n)
rprim2[ :,2] = rprim0[ :,2] * (1.00 + eps * k)

namei = int(round(1000 * (1.00 + eps * i)))
namej = int(round(1000 * (1.00 + eps * j)))
namek = int(round(1000 * (1.00 + eps * k)))
namel = int(round(1000 * (1.00 + eps * l)))
namem = int(round(1000 * (1.00 + eps * m)))
namen = int(round(1000 * (1.00 + eps * n)))
formatted_namei = f"{namei:04d}_{namej:04d}_{namek:04d}_{namel:04d}_{namem:04d}_{namen:04d}"

structure2=structure.copy()
structure2.lattice=Lattice(rprim2)
structures_new[formatted_namei] = structure2

return structures_new
elif 3 <= spgrp_number <= 15:
disp=[[1,1,1,1], [0,1,1,1], [2,1,1,1], [1,0,1,1], [1,2,1,1], [1,1,0,1], [1,1,2,1], [1,1,1,0],
[1,1,1,2], [0,0,1,1], [1,0,0,1], [1,1,0,0], [0,1,0,1], [1,0,1,0], [0,1,1,0]]
if abs(rprim[1, 0]) > 1e-9 or abs(rprim[0, 1]) > 1e-9 or abs(rprim[2, 1]) > 1e-9 or abs(rprim[1, 2]) > 1e-9:
print("Error: Monoclinic structure with yx=xy=0 and yz=zy=0 lattice required.")
elif abs(rprim[0, 2]) > 1e-9 :
print("Warning: The lattice is oriented such that xz = 0.")
rprim0 = np.copy(rprim)
a=rprim[0, :]
b=rprim[1, :]
c=rprim[2, :]
norm_a = np.linalg.norm(a)
norm_b = np.linalg.norm(b)
norm_c = np.linalg.norm(c)

# Compute angles between vectors
cos_ab = np.dot(a, b) / (norm_a * norm_b)
cos_ac = np.dot(a, c) / (norm_a * norm_c)
cos_bc = np.dot(b, c) / (norm_b * norm_c)

rprim0[0,0] = norm_a
rprim0[0,2] = 0.0
rprim0[1,1] = norm_b
rprim0[2,0] = norm_c*cos_ac
rprim0[2,2] = norm_c*np.sqrt(1-cos_ac**2)
print("Old rprim:")
print(rprim)
print("New rprim:")
print(rprim0)

for pair in disp:
i,j,k,l = pair
rprim2[ :,0] = rprim0[ :,0] * (1.00 + eps * i) +rprim0[ :,2] * (eps * l)
rprim2[ :,1] = rprim0[ :,1] * (1.00 + eps * j)
rprim2[ :,2] = rprim0[ :,2] * (1.00 + eps * k)

namei = int(round(1000 * (1.00 + eps * i)))
namej = int(round(1000 * (1.00 + eps * j)))
namek = int(round(1000 * (1.00 + eps * k)))
namel = int(round(1000 * (1.00 + eps * l)))
formatted_namei = f"{namei:04d}_{namej:04d}_{namek:04d}_{namel:04d}"

structure2=structure.copy()
structure2.lattice=Lattice(rprim2)
structures_new[formatted_namei] = structure2

return structures_new
elif 16 <= spgrp_number <= 74:
disp=[[0,0,1],[0,1,0],[1,0,0],[1,1,1],[0,1,1],[2,1,1],[1,0,1],[1,2,1],[1,1,0],[1,1,2]]
for pair in disp:
i,j,k = pair
rprim2[ :,0] = rprim[ :,0] * (1.00 + eps * i)
rprim2[ :,1] = rprim[ :,1] * (1.00 + eps * j)
rprim2[ :,2] = rprim[ :,2] * (1.00 + eps * k)

namei = int(round(1000 * (1.00 + eps * i)))
namej = int(round(1000 * (1.00 + eps * j)))
namek = int(round(1000 * (1.00 + eps * k)))
formatted_namei = f"{namei:04d}_{namej:04d}_{namek:04d}"

structure2=structure.copy()
structure2.lattice=Lattice(rprim2)
structures_new[formatted_namei] = structure2

return structures_new
elif 75 <= spgrp_number <= 194:
disp=[[0,0],[1,1],[0,1],[2,1],[1,0],[1,2]]
for pair in disp:
i, k = pair
rprim2[ :,0] = rprim[ :,0] * (1.00 + eps * i)
rprim2[ :,1] = rprim[ :,1] * (1.00 + eps * i)
rprim2[ :,2] = rprim[ :,2] * (1.00 + eps * k)

namei = int(round(1000 * (1.00 + eps * i)))
namek = int(round(1000 * (1.00 + eps * k)))
formatted_namei = f"{namei:04d}_{namek:04d}"
rprim_new[formatted_namei] = rprim2

structure2=structure.copy()
structure2.lattice=Lattice(rprim2)
structures_new[formatted_namei] = structure2

return structures_new
elif 195 <= spgrp_number <= 230:
for i in range(3):
rprim2[ :,0] = rprim[ :,0] * (1.00 + eps * i)
rprim2[ :,1] = rprim[ :,1] * (1.00 + eps * i)
rprim2[ :,2] = rprim[ :,2] * (1.00 + eps * i)
namei = int(round(1000 * (1.00 + eps * i)))
formatted_namei = f"{namei:04d}"

structure2=structure.copy()
structure2.lattice=Lattice(rprim2)
structures_new[formatted_namei] = structure2
return structures_new

Loading

0 comments on commit 567fa3f

Please sign in to comment.