Python Molecule.xyzs方法代码示例

# 需要导入模块: from molecule import Molecule [as 别名]
# 或者: from molecule.Molecule import xyzs [as 别名]
    def fgrad(x, indicate = False):
        """ Calculate the objective function and its derivatives. """
        # If the optimization algorithm tries to calculate twice for the same point, do nothing.
        # if x == fgrad.x0: return

        xyz = independent_vars_to_xyz(x)
        # these methods require 3d input
        xyzlist = np.array([xyz])
        my_bonds = core.bonds(xyzlist, ibonds).flatten()
        my_angles = core.angles(xyzlist, iangles).flatten()
        my_dihedrals = core.dihedrals(xyzlist, idihedrals, anchor=dihedrals).flatten()

        # Deviations of internal coordinates from ideal values.
        d1 = w1*(my_bonds - bonds)
        d2 = w2*(my_angles - angles)
        d3 = w3*(my_dihedrals - dihedrals)

        # Include an optional term if we have an anchor point.
        if xrefi != None:
            d4 = (x - xrefi).flatten() * w1 * w_xref
            fgrad.error = np.r_[d1, d2, np.arctan2(np.sin(d3), np.cos(d3)), d4]
        else:
            fgrad.error = np.r_[d1, d2, np.arctan2(np.sin(d3), np.cos(d3))]

        # The objective function contains another contribution from the Morse potential.
        fgrad.X = np.dot(fgrad.error, fgrad.error)
        d1s = np.dot(d1, d1)
        d2s = np.dot(d2, d2)
        d3s = np.dot(d3, d3)
        M = Molecule()
        M.elem = elem
        M.xyzs = [np.array(xyz)*10]
        if w_morse != 0.0:
            EMorse, GMorse = PairwiseMorse(M)
            EMorse = EMorse[0]
            GMorse = GMorse[0]
        else:
            EMorse = 0.0
            GMorse = np.zeros((n_atoms, 3), dtype=float)
        if indicate: 
            if fgrad.X0 != None:
                print ("LSq: %.4f (%+.4f) Distance: %.4f (%+.4f) Angle: %.4f (%+.4f) Dihedral: %.4f (%+.4f) Morse: % .4f (%+.4f)" % 
                       (fgrad.X, fgrad.X - fgrad.X0, d1s, d1s - fgrad.d1s0, d2s, d2s - fgrad.d2s0, d3s, d3s - fgrad.d3s0, EMorse, EMorse - fgrad.EM0)), 
            else:
                print "LSq: %.4f Distance: %.4f Angle: %.4f Dihedral: %.4f Morse: % .4f" % (fgrad.X, d1s, d2s, d3s, EMorse), 
                fgrad.X0 = fgrad.X
                fgrad.d1s0 = d1s
                fgrad.d2s0 = d2s
                fgrad.d3s0 = d3s
                fgrad.EM0 = EMorse
        fgrad.X += w_morse*EMorse

        # Derivatives of internal coordinates w/r.t. Cartesian coordinates.
        d_bonds = core.bond_derivs(xyz, ibonds) * w1
        d_angles = core.angle_derivs(xyz, iangles) * w2
        d_dihedrals = core.dihedral_derivs(xyz, idihedrals) * w3

        if xrefi != None:
            # the derivatives of the internal coordinates wrt the cartesian
            # this is 2d, with shape equal to n_internal x n_cartesian
            d_internal = np.vstack([gxyz_to_independent_vars(d_bonds.reshape((len(ibonds), -1))),
                                    gxyz_to_independent_vars(d_angles.reshape((len(iangles), -1))),
                                    gxyz_to_independent_vars(d_dihedrals.reshape((len(idihedrals), -1))),
                                    np.eye(len(x)) * w1 * w_xref])
        else:
            # the derivatives of the internal coordinates wrt the cartesian
            # this is 2d, with shape equal to n_internal x n_cartesian
            d_internal = np.vstack([gxyz_to_independent_vars(d_bonds.reshape((len(ibonds), -1))),
                                    gxyz_to_independent_vars(d_angles.reshape((len(iangles), -1))),
                                    gxyz_to_independent_vars(d_dihedrals.reshape((len(idihedrals), -1)))])

        # print fgrad.error.shape, d_internal.shape
        # print d_internal.shape
        # print xyz_to_independent_vars(d_internal).shape
        fgrad.G = 2*np.dot(fgrad.error, d_internal)
        fgrad.G += xyz_to_independent_vars(w_morse*GMorse.flatten())