Python units.Bohr方法代碼示例

# 需要導入模塊: from ase import units [as 別名]
# 或者: from ase.units import Bohr [as 別名]
def get_hessian(atoms, delta=0.005, mass_weighted=False):
    """
    Calculate (mass weighted) hessian using central diff formula.
    :param atoms: atoms object with defined calculator
    :param delta: step size for numeric differentiation
    :type atoms: ase.Atoms
    :type delta: float
    :return: numpy square symmetric array
    """
    # convert delta to Angs
    delta *= Bohr
    # allocate matrix
    l = len(atoms)
    H = np.zeros((3 * l, 3 * l), dtype=np.float64)
    r = 0
    # gradients matrix
    for i, j in product(range(l), range(3)):
        g = np.zeros((l, 3))
        for k in (-1, 1):
            atoms1 = atoms.copy()
            atoms1[i].position[j] += k * delta
            atoms1.set_calculator(atoms.get_calculator())
            g += - k * atoms1.get_forces()

        H[r] = 0.5 * g.flatten()
        r += 1
    # check symmetry assuming gradients computed with 10^-3 Hartree/Bohr precision
    gprec = 0.001 * Hartree
    assert np.max(np.abs(H - H.T)) < gprec, np.max(np.abs(H - H.T))
    # Hessian
    H /= delta
    # symmetrize
    H = 0.5 * (H + H.T)
    # mass weight
    if mass_weighted:
        v = np.sqrt(atoms.get_masses()).repeat(3).reshape(-1, 1)
        H /= np.dot(v, v.T)
    return H 

# 需要導入模塊: from ase import units [as 別名]
# 或者: from ase.units import Bohr [as 別名]
def format_geo(self, is_angstrom=False):
    '''Converts geo_info and geo_spec to a universal format.
    **Parameters:**

    is_angstrom : bool, optional
      If True, input is assumed to be in Angstrom and positions are converted to Bohr radii.
    '''
    for i in self.geo_info:
      i[0] = get_atom_symbol(i[0])
      i[2] = float(i[-1])
    self.geo_info = numpy.array(self.geo_info)
    self.geo_spec = numpy.array(self.geo_spec,dtype=float)
    if is_angstrom:
      self.geo_spec *= aa_to_a0 

# 需要導入模塊: from ase import units [as 別名]
# 或者: from ase.units import Bohr [as 別名]
def get_ase_atoms(self,bbox=None,**kwargs):
    '''Create an ASE atoms object.
    (cf. https://wiki.fysik.dtu.dk/ase/ase/atoms.html )

    **Parameters:**

    bbox : list of floats (bbox=[xmin,xmax,ymin,ymax,zmin,zmax]), optional
      If not None, sets the unit cell to the grid boundaries and moves the
      molecule in its center.

    **Returns:**

    atoms : Atoms object
      See https://wiki.fysik.dtu.dk/ase/ase/atoms.html for details

    .. Note::

      ASE has to be in the PYTHONPATH
    '''
    from ase import Atoms
    from ase.units import Bohr

    atoms = Atoms("".join(self.geo_info[:,0]),
                  positions=self.geo_spec*Bohr,
                  **kwargs)
    if bbox is not None:
      if len(bbox) != 6:
        raise ValueError("bbox has to have 6 elements")
      bbox = numpy.array(bbox)
      atoms.translate(-bbox[::2]*Bohr)
      atoms.cell = numpy.eye(3) * (bbox[1::2] - bbox[::2])*Bohr

    return atoms
  # Synonym 

# 需要導入模塊: from ase import units [as 別名]
# 或者: from ase.units import Bohr [as 別名]
def load_data(dbpath):
    logging.info('Downloading GDB-9 data...')
    tmpdir = tempfile.mkdtemp('gdb9')
    tar_path = os.path.join(tmpdir, 'gdb9.tar.gz')
    raw_path = os.path.join(tmpdir, 'gdb9_xyz')
    url = 'https://ndownloader.figshare.com/files/3195389'

    try:
        urllib.request.urlretrieve(url, tar_path)
        logging.info("Done.")
    except HTTPError as e:
        logging.error("HTTP Error:", e.code, url)
        return False
    except URLError as e:
        logging.error("URL Error:", e.reason, url)
        return False

    tar = tarfile.open(tar_path)
    tar.extractall(raw_path)
    tar.close()

    prop_names = ['rcA', 'rcB', 'rcC', 'mu', 'alpha', 'homo', 'lumo',
                  'gap', 'r2', 'zpve', 'energy_U0', 'energy_U', 'enthalpy_H',
                  'free_G', 'Cv']
    conversions = [1., 1., 1., 1., Bohr ** 3 / Ang ** 3,
                   Hartree / eV, Hartree / eV, Hartree / eV,
                   Bohr ** 2 / Ang ** 2, Hartree / eV,
                   Hartree / eV, Hartree / eV, Hartree / eV,
                   Hartree / eV, 1.]

    logging.info('Parse xyz files...')
    with connect(dbpath) as con:
        for i, xyzfile in enumerate(os.listdir(raw_path)):
            xyzfile = os.path.join(raw_path, xyzfile)

            if i % 10000 == 0:
                logging.info('Parsed: ' + str(i) + ' / 133885')
            properties = {}
            tmp = os.path.join(tmpdir, 'tmp.xyz')

            with open(xyzfile, 'r') as f:
                lines = f.readlines()
                l = lines[1].split()[2:]
                for pn, p, c in zip(prop_names, l, conversions):
                    properties[pn] = float(p) * c
                with open(tmp, "wt") as fout:
                    for line in lines:
                        fout.write(line.replace('*^', 'e'))

            with open(tmp, 'r') as f:
                ats = list(read_xyz(f, 0))[0]

            con.write(ats, key_value_pairs=properties)
    logging.info('Done.')

    return True