Source code for dpgen.auto_test.lib.siesta

import numpy as np
from dpdata.periodic_table import Element


def _make_siesta_01_common(sys_data, ecut, ediff, mixingWeight, NumberPulay):
    tot_natoms = sum(sys_data["atom_numbs"])
    ntypes = len(sys_data["atom_names"])
    ret = ""
    ret += "SystemName        system\n"
    ret += "SystemLabel       system\n"
    ret += "NumberOfAtoms     %d" % tot_natoms
    ret += "\nNumberOfSpecies   %d\n" % ntypes
    ret += "\n"
    ret += "WriteForces       T\n"
    ret += "WriteCoorStep     T\n"
    ret += "WriteCoorXmol     T\n"
    ret += "WriteMDXmol       T\n"
    ret += "WriteMDHistory    T\n\n"

    ret += "MeshCutoff            %s" % str(ecut)
    ret += " Ry\n"
    ret += "DM.MixingWeight       %f\n" % mixingWeight
    ret += "DM.Tolerance          %e\n" % ediff
    ret += "DM.UseSaveDM          true\n"
    ret += "DM.NumberPulay         %d\n" % NumberPulay
    ret += "MD.UseSaveXV           T\n\n"

    ret += "XC.functional          GGA\n"
    ret += "XC.authors             PBE\n"
    ret += "DM.UseSaveDM           F\n"
    ret += "WriteDM                F\n"
    ret += "WriteDM.NetCDF         F\n"
    ret += "WriteDMHS.NetCDF       F\n"
    return ret


def _make_siesta_02_species(sys_data, pps):
    atom_nums = sys_data["atom_numbs"]
    atom_names = sys_data["atom_names"]
    ntypes = len(atom_nums)
    assert ntypes == len(atom_names)
    assert ntypes == len(pps)
    ret = ""
    ret += "%block Chemical_Species_label\n"
    for i in range(0, len(atom_names)):
        ret += (
            str(i + 1)
            + "\t"
            + str(Element(atom_names[i]).Z)
            + "\t"
            + atom_names[i]
            + "\n"
        )
    ret += "%endblock Chemical_Species_label\n"
    return ret


# ## kpoints !!! can not understand
def _make_siesta_03_kpoint(sys_data, kspacing):
    cell = sys_data["cells"][0]
    cell = np.reshape(cell, [3, 3])
    ## np.linalg.inv()：矩阵求逆
    rcell = np.linalg.inv(cell)
    ## .T 矩阵转置
    rcell = rcell.T
    # np.ceil()是向上取整，与四舍五入无关 -5.6 --> -5
    # np.linalg.norm：进行范数运算，范数是对向量（或者矩阵）的度量，是一个标量（scalar）
    kpoints = [
        (np.ceil(2 * np.pi * np.linalg.norm(ii) / kspacing).astype(int)) for ii in rcell
    ]
    ret = ""
    ret += "%block kgrid_Monkhorst_Pack\n"
    ret += "%d" % kpoints[0]
    ret += "\t0\t0\t0.0\n"

    ret += "0\t"
    ret += "%d" % kpoints[1]
    ret += "\t0\t0.0\n"

    ret += "0\t0\t"
    ret += "%d" % kpoints[2]
    ret += "\t0.0\n"

    ret += "%endblock kgrid_Monkhorst_Pack\n"
    return ret


### coordinate
def _make_siesta_04_ucVectorCoord(sys_data):
    cell = sys_data["cells"][0]
    cell = np.reshape(cell, [3, 3])
    coordinates = sys_data["coords"][0]
    atom_names = sys_data["atom_names"]
    atom_numbs = sys_data["atom_numbs"]
    ntypes = len(atom_names)
    ret = ""
    ret += "LatticeConstant try_input_output.00 Ang\n"
    ret += "%block LatticeVectors\n"
    for ii in range(3):
        for jj in range(3):
            ret += "%f " % cell[ii][jj]
        ret += "\n"
    ret += "%endblock LatticeVectors\n"

    ret += "\n"
    ret += "AtomicCoordinatesFormat Ang\n"
    ret += "%block AtomicCoordinatesAndAtomicSpecies\n"
    cc = 0
    for ii in range(ntypes):
        for jj in range(atom_numbs[ii]):
            ret += "%f %f %f %d %s\n" % (
                coordinates[cc][0],
                coordinates[cc][1],
                coordinates[cc][2],
                ii + 1,
                atom_names[ii],
            )
            cc += 1
    ret += "%endblock AtomicCoordinatesAndAtomicSpecies"
    return ret




def make_siesta_input(sys_data, fp_pp_files, fp_params):
    ecut = fp_params["ecut"]
    ediff = fp_params["ediff"]
    mixingWeight = fp_params["mixingWeight"]
    NumberPulay = fp_params["NumberPulay"]
    kspacing = fp_params["kspacing"]
    ret = ""
    ret += _make_siesta_01_common(sys_data, ecut, ediff, mixingWeight, NumberPulay)
    ret += "\n"
    ret += _make_siesta_02_species(sys_data, fp_pp_files)
    ret += "\n"
    ret += _make_siesta_03_kpoint(sys_data, kspacing)
    ret += "\n"
    ret += _make_siesta_04_ucVectorCoord(sys_data)
    ret += "\n"
    return ret