upload elastic

_book/examples/elastic/INPUT

@@ -0,0 +1,24 @@
+INPUT_PARAMETERS
+#Parameters (1.General)
+suffix          C8
+calculation     relax
+esolver_type    ksdft
+symmetry        1
+pseudo_dir      ../
+orbital_dir     ../
+cal_stress      1
+cal_force       1
+#Parameters (2.Iteration)
+ecutwfc         100
+scf_thr         1e-7
+scf_nmax        50
+#Parameters (3.Basis)
+basis_type      lcao
+# kpar            8
+gamma_only      0
+#Parameters (4.Smearing)
+smearing_method gaussian
+smearing_sigma  0.002
+#Parameters (5.Mixing)
+mixing_type     broyden
+mixing_beta     0.7

_book/examples/elastic/KPT

@@ -0,0 +1,4 @@
+K_POINTS
+0
+MP
+14 14 14 0 0 0

_book/examples/elastic/compute_dfm.py

@@ -0,0 +1,110 @@
+from monty.serialization import loadfn, dumpfn
+from pymatgen.analysis.elasticity.elastic import ElasticTensor
+from pymatgen.analysis.elasticity.stress import Stress
+import glob, os, sys
+import numpy as np
+
+def get_stress_vasp(lines: str) -> np.ndarray[3, 3]:
+    stress = np.zeros([3,3])
+    for line in lines:
+        if "in kB" in line:
+            stress_xx = float(line.split()[2])
+            stress_yy = float(line.split()[3])
+            stress_zz = float(line.split()[4])
+            stress_xy = float(line.split()[5])
+            stress_yz = float(line.split()[6])
+            stress_zx = float(line.split()[7])
+            stress[0] = [stress_xx, stress_xy, stress_zx]
+            stress[1] = [stress_xy, stress_yy, stress_yz]
+            stress[2] = [stress_zx, stress_yz, stress_zz]
+    return stress
+
+def get_stress_abacus(lines: str) -> np.ndarray[3, 3]:
+    stress = np.zeros([3,3])
+    for idx, line in enumerate(lines):
+        if "TOTAL-STRESS (KBAR)" in line:
+            stress_xx = float(lines[idx+2].split()[0])
+            stress_yy = float(lines[idx+3].split()[1])
+            stress_zz = float(lines[idx+4].split()[2])
+            stress_xy = float(lines[idx+2].split()[1])
+            stress_yz = float(lines[idx+3].split()[2])
+            stress_zx = float(lines[idx+2].split()[2])
+            stress[0] = [stress_xx, stress_xy, stress_zx]
+            stress[1] = [stress_xy, stress_yy, stress_yz]
+            stress[2] = [stress_zx, stress_yz, stress_zz]
+    return stress
+
+try:
+    run_type = sys.argv[1]
+except:
+    print("Usage: python compute_dfm.py [abacus|vasp]")
+    sys.exit(1)
+
+if  run_type == "abacus":
+    OUTCAR = "OUT.*/running_*.log"
+elif run_type == "vasp":
+    OUTCAR = "OUTCAR"
+
+cwd = os.getcwd()
+# print(cwd)
+
+# equi_stress
+equi = glob.glob(os.path.join(cwd, "relax/", OUTCAR))[0]
+# print(equi)
+with open(equi, "r") as fin:
+    lines = fin.read().split("\n")
+if run_type == "abacus":
+    equi_stress = Stress(get_stress_abacus(lines) * (-1000))
+elif run_type == "vasp":
+    equi_stress = Stress(get_stress_vasp(lines) * (-1000))
+    # print(equi_stress)
+
+
+# read all the task dir
+task_dirs = glob.glob("task.*")
+lst_strain = []
+lst_stress = []
+for ii in task_dirs:
+    os.chdir(os.path.join('./', ii))
+
+    strain = loadfn("strain.json")
+    # print(strain, strain.shape)
+
+    stress = np.zeros([3,3])
+    OUTCAR = glob.glob(OUTCAR)[0]
+    with open(OUTCAR, "r") as fin:
+        lines = fin.read().split("\n")
+    if run_type == "abacus":
+        stress = get_stress_abacus(lines)
+    elif run_type == "vasp":
+        stress = get_stress_vasp(lines)
+    # print(stress, stress.shape)
+    os.chdir(cwd)
+    lst_strain.append(strain)
+    lst_stress.append(Stress(stress * (-1000)))
+
+# print(lst_strain)
+et = ElasticTensor.from_independent_strains(lst_strain, lst_stress, eq_stress=equi_stress, vasp=False)
+
+res_data = {}
+ptr_data = '# Elastic Constants in GPa\n'
+res_data["elastic_tensor"] = []
+for ii in range(6):
+    for jj in range(6):
+        res_data["elastic_tensor"].append(et.voigt[ii][jj] / 1e4)
+        ptr_data += "%7.2f " % (et.voigt[ii][jj] / 1e4)
+    ptr_data += "\n"
+BV = et.k_voigt / 1e4
+GV = et.g_voigt / 1e4
+EV = 9 * BV * GV / (3 * BV + GV)
+uV = 0.5 * (3 * BV - 2 * GV) / (3 * BV + GV)
+res_data["BV"] = BV
+res_data["GV"] = GV
+res_data["EV"] = EV
+res_data["uV"] = uV
+ptr_data += "# Bulk   Modulus BV = %.2f GPa\n" % BV
+ptr_data += "# Shear  Modulus GV = %.2f GPa\n" % GV
+ptr_data += "# Youngs Modulus EV = %.2f GPa\n" % EV
+ptr_data += "# Poission Ratio uV = %.2f " % uV
+print(ptr_data)
+dumpfn(res_data, "elastic.json", indent=4)

_book/examples/elastic/gene_dfm.py

@@ -0,0 +1,96 @@
+from pymatgen.core.structure import Structure
+from pymatgen.analysis.elasticity.elastic import Strain
+from pymatgen.analysis.elasticity.strain import DeformedStructureSet
+import os, sys, dpdata, glob
+from monty.serialization import dumpfn
+
+try:
+    run_type = sys.argv[1]
+except:
+    print("Usage: python gene_dfm.py [abacus|vasp]")
+    sys.exit(1)
+
+cwd = os.getcwd()
+path_to_equi = os.path.join(cwd, 'relax')
+
+if  run_type == "abacus":
+    CONTCAR = os.path.join('OUT.*', 'STRU_ION_D')
+    POSCAR = "STRU"
+    INCAR = "INPUT"
+    KPOINTS = "KPT"
+elif run_type == "vasp":
+    CONTCAR = "CONTCAR"
+    POSCAR = "POSCAR"
+    INCAR = "INCAR"
+    KPOINTS = "KPOINTS"
+
+# print(CONTCAR)
+equi_contcar = glob.glob(os.path.join(path_to_equi, CONTCAR))[0]
+# print(equi_contcar)
+if not os.path.exists(equi_contcar):
+    raise RuntimeError("Please do relaxation first!")
+
+if run_type == "abacus":
+    stru = dpdata.System(equi_contcar, fmt = "stru")
+    stru.to("poscar", "POSCAR.tmp")
+    ss = Structure.from_file("POSCAR.tmp")
+    os.remove("POSCAR.tmp")
+elif run_type == "vasp":
+    ss = Structure.from_file(equi_contcar)
+
+norm_strains = [-0.010, -0.005, 0.005, 0.010]
+shear_strains = [-0.010, -0.005, 0.005, 0.010]
+
+dfm_ss = DeformedStructureSet(ss, symmetry=False, norm_strains=norm_strains, shear_strains=shear_strains)
+# print(dfm_ss)
+n_dfm = len(dfm_ss)
+
+print("gen with norm " + str(norm_strains))
+print("gen with shear " + str(shear_strains))
+for ii in range(n_dfm):
+    output_task = os.path.join('./', "task.%03d" % ii)
+    os.makedirs(output_task, exist_ok=True)
+    os.chdir(output_task)
+    dfm_ss.deformed_structures[ii].to("POSCAR", fmt = "POSCAR")
+    if run_type == "abacus":
+        stru = dpdata.System("POSCAR", fmt="vasp/poscar")
+        n_atoms = len(stru["atom_names"])
+        atom_mass = []
+        pseudo = []
+        orb = []
+        with open(equi_contcar, "r") as f:
+            lines = f.readlines()
+        for idx, line in enumerate(lines):
+            if "ATOMIC_SPECIES" in line:
+                for i in range(n_atoms):
+                    atom_mass.append(float(lines[idx+i+1].split()[1]))
+                    pseudo.append(lines[idx+i+1].split()[2])
+            if "NUMERICAL_ORBITAL" in line:
+                for i in range(n_atoms):
+                    orb.append(lines[idx+i+1])
+        if orb == []:
+            stru.to("stru", "STRU", mass=atom_mass, pp_file=pseudo)
+        else:
+            stru.to("stru", "STRU", mass=atom_mass, pp_file=pseudo, numerical_orbital=orb)
+        os.remove("POSCAR")
+    os.system("cp ../{} .".format(INCAR))
+    if run_type == "abacus":
+        with open(INCAR, "r") as f:
+            lines = f.readlines()
+        pseudo_dir = "../"
+        orb_dir = "../"
+        for line in lines:
+            if 'pseudo_dir' in line:
+                if pseudo_dir != line.split()[1].strip():
+                    line = line.replace(line, 'pseudo_dir ' + pseudo_dir)
+            if 'orb_dir' in line:
+                if orb_dir != line.split()[1].strip():
+                    line = line.replace(line, 'orb_dir ' + orb_dir)
+        with open(INCAR, "w") as f:
+            f.writelines(lines)
+    os.system("cp ../{} .".format(KPOINTS))
+    if run_type == "vasp":
+        os.system("cp ../POTCAR .")
+    df = Strain.from_deformation(dfm_ss.deformations[ii])
+    dumpfn(df.as_dict(), "strain.json", indent=4)
+    os.chdir(cwd)

_book/examples/elastic/relax/INPUT

@@ -0,0 +1,24 @@
+INPUT_PARAMETERS
+#Parameters (1.General)
+suffix          C8
+calculation     relax
+esolver_type    ksdft
+symmetry        1
+pseudo_dir      ../
+orbital_dir     ../
+cal_stress      1
+cal_force       1
+#Parameters (2.Iteration)
+ecutwfc         100
+scf_thr         1e-7
+scf_nmax        50
+#Parameters (3.Basis)
+basis_type      lcao
+# kpar            8
+gamma_only      0
+#Parameters (4.Smearing)
+smearing_method gaussian
+smearing_sigma  0.002
+#Parameters (5.Mixing)
+mixing_type     broyden
+mixing_beta     0.7

_book/examples/elastic/relax/KPT

@@ -0,0 +1,4 @@
+K_POINTS
+0
+MP
+14 14 14 0 0 0

_book/examples/elastic/relax/STRU

@@ -0,0 +1,29 @@
+ATOMIC_SPECIES
+C 12.011 C_ONCV_PBE-1.0.upf upf201
+
+NUMERICAL_ORBITAL
+C_gga_7au_100Ry_2s2p1d.orb
+
+LATTICE_CONSTANT
+1.889726
+
+LATTICE_VECTORS
+   3.5736050709821550    0.0000000000000000    0.0000000000000000
+   0.0000000000000000    3.5736050709821550    0.0000000000000000
+   0.0000000000000000    0.0000000000000000    3.5736050709821550
+
+ATOMIC_POSITIONS
+Direct
+
+C #label
+0 #magnetism
+8 #number of atoms
+   0.2500000000000000    0.2500000000000000    0.2500000000000000
+   0.0000000000000000    0.0000000000000000    0.0000000000000000
+   0.2500000000000000    0.7500000000000000    0.7500000000000000
+   0.0000000000000000    0.5000000000000000    0.5000000000000000
+   0.7500000000000000    0.2500000000000000    0.7500000000000000
+   0.5000000000000000    0.0000000000000000    0.5000000000000000
+   0.7500000000000000    0.7500000000000000    0.2500000000000000
+   0.5000000000000000    0.5000000000000000    0.0000000000000000
+

_book/examples/elastic/requirements.txt

@@ -0,0 +1,4 @@
+dpdata
+monty
+numpy
+pymatgen

_book/examples/elastic/run_task.sh

@@ -0,0 +1,11 @@
+#!/bin/bash
+for i in task.*
+do
+cd ./$i
+pwd
+# 超算环境注意修改sub.sh中的内容
+sbatch ../sub.sh
+# 或者直接运行abacus
+# mpirun -n 4 abacus
+cd ../
+done

_book/examples/elastic/sub.sh

@@ -0,0 +1,7 @@
+#!/bin/bash
+#SBATCH -p master
+#SBATCH -J lcao-elas
+#SBATCH -n 6
+
+intel
+mpirun -n 6 abacus | tee abacus.out

abacus-upf.md

@@ -412,7 +412,7 @@ upfconv.x -u  al.upf #利用QE的upflib将UPFv.1转化成UPFv.2格式
 
 ## 3. ld1.x
 
-QE 的 atomic 模块中的 `ld1.x` 支持生成赝势。其不仅可以生成模守恒赝势、还支持超软赝势、PAW 方法，支持全相对论、标量（非）相对论赝势（rel, non-rel/sca-rel），其径向波函数支持 [TM](https://journals.aps.org/prb/abstract/10.1103/PhysRevB.43.1993)（更稳定）与 [RRKJ](https://journals.aps.org/prb/abstract/10.1103/PhysRevB.41.1227) 两种方法，支持交换关联近似 7 类 LDA（ pz）,GGA（pbe, pbesol, revpbe, bp, wc, pw91。
+QE 的 atomic 模块中的 `ld1.x` 支持生成赝势。其不仅可以生成模守恒赝势、还支持超软赝势、PAW 方法，支持全相对论、标量（非）相对论赝势（rel, non-rel/sca-rel），其径向波函数支持 [TM](https://journals.aps.org/prb/abstract/10.1103/PhysRevB.43.1993)（更稳定）与 [RRKJ](https://journals.aps.org/prb/abstract/10.1103/PhysRevB.41.1227) 两种方法，支持交换关联近似 7 类 LDA（pz）,GGA（pbe, pbesol, revpbe, bp, wc, pw91）。
 
 下载好 QE 软件后，可以直接安装：
 
@@ -497,7 +497,7 @@ element='C Si Ge'
 /
 &inputp
    pseudotype=2,
-   file_pseudopw='Al.pbe-n-nc.UPF', #输入赝势
+   file_pseudopw='Al.pbe-n-nc.UPF', #输出赝势
    lloc=2,
    nlcc=.true.,
    tm=.true.