How to use the forcebalance.finite_difference.f12d3p function in forcebalance

    @param[in] pdb OpenMM PDB object
    @param[in] FF ForceBalance force field object
    @param[in] xyzs List of OpenMM positions
    @param[in] settings OpenMM settings for creating the System
    @param[in] boxes Periodic box vectors
    @return G First derivative of the energies in a N_param x N_coord array

    """

    G        = np.zeros((FF.np,len(xyzs)))
    if not AGrad:
        return G
    E0       = energy_driver(mvals, pdb, FF, xyzs, settings, simulation, boxes)
    CheckFDPts = False
    for i in range(FF.np):
        G[i,:], _ = f12d3p(fdwrap(energy_driver,mvals,i,key=None,pdb=pdb,FF=FF,xyzs=xyzs,settings=settings,simulation=simulation,boxes=boxes),h,f0=E0)
        if CheckFDPts:
            # Check whether the number of finite difference points is sufficient.  Forward difference still gives residual error of a few percent.
            G1 = f1d7p(fdwrap(energy_driver,mvals,i,key=None,pdb=pdb,FF=FF,xyzs=xyzs,settings=settings,simulation=simulation,boxes=boxes),h)
            dG = G1 - G[i,:]
            dGfrac = (G1 - G[i,:]) / G[i,:]
            print "Parameter %3i 7-pt vs. central derivative : RMS, Max error (fractional) = % .4e % .4e (% .4e % .4e)" % (i, np.sqrt(np.mean(dG**2)), max(np.abs(dG)), np.sqrt(np.mean(dGfrac**2)), max(np.abs(dGfrac)))
    return G

   @todo This is a sufficiently general function to be merged into openmmio.py?
   @param[in] mvals Mathematical parameter values
   @param[in] pdb OpenMM PDB object
   @param[in] FF ForceBalance force field object
   @param[in] xyzs List of OpenMM positions
   @return G First derivative of the energies in a N_param x N_coord array
   @return Hd Second derivative of the energies (i.e. diagonal Hessian elements) in a N_param x N_coord array

   """

   G        = np.zeros((FF.np,len(xyzs)))
   Hd       = np.zeros((FF.np,len(xyzs)))
   E0       = energy_driver(mvals, pdb, FF, xyzs, boxes)
   for i in range(FF.np):
      G[i,:], Hd[i,:] = f12d3p(fdwrap(energy_driver,mvals,i,key=None,pdb=pdb,FF=FF,xyzs=xyzs,boxes=boxes),h,f0=E0)
   return G, Hd

    @param[in] settings OpenMM settings for creating the System
    @param[in] boxes Periodic box vectors
    @return G First derivative of the energies in a N_param x N_coord array

    """

    G        = np.zeros((FF.np,len(xyzs)))
    GDx      = np.zeros((FF.np,len(xyzs)))
    GDy      = np.zeros((FF.np,len(xyzs)))
    GDz      = np.zeros((FF.np,len(xyzs)))
    if not AGrad:
        return G, GDx, GDy, GDz
    ED0      = energy_driver(mvals, pdb, FF, xyzs, settings, simulation, boxes, dipole=True)
    CheckFDPts = False
    for i in range(FF.np):
        EDG, _   = f12d3p(fdwrap(energy_driver,mvals,i,key=None,pdb=pdb,FF=FF,xyzs=xyzs,settings=settings,simulation=simulation,boxes=boxes,dipole=True),h,f0=ED0)
        G[i,:]   = EDG[:,0]
        GDx[i,:] = EDG[:,1]
        GDy[i,:] = EDG[:,2]
        GDz[i,:] = EDG[:,3]
    return G, GDx, GDy, GDz

    @param[in] pdb OpenMM PDB object
    @param[in] FF ForceBalance force field object
    @param[in] xyzs List of OpenMM positions
    @param[in] settings OpenMM settings for creating the System
    @param[in] boxes Periodic box vectors
    @return G First derivative of the energies in a N_param x N_coord array

    """

    G        = np.zeros((FF.np,len(xyzs)))
    if not AGrad:
        return G
    E0       = energy_driver(mvals, pdb, FF, xyzs, settings, simulation, boxes)
    CheckFDPts = False
    for i in range(FF.np):
        G[i,:], _ = f12d3p(fdwrap(energy_driver,mvals,i,key=None,pdb=pdb,FF=FF,xyzs=xyzs,settings=settings,simulation=simulation,boxes=boxes),h,f0=E0)
        if CheckFDPts:
            # Check whether the number of finite difference points is sufficient.  Forward difference still gives residual error of a few percent.
            G1 = f1d7p(fdwrap(energy_driver,mvals,i,key=None,pdb=pdb,FF=FF,xyzs=xyzs,settings=settings,simulation=simulation,boxes=boxes),h)
            dG = G1 - G[i,:]
            dGfrac = (G1 - G[i,:]) / G[i,:]
            print "Parameter %3i 7-pt vs. central derivative : RMS, Max error (fractional) = % .4e % .4e (% .4e % .4e)" % (i, np.sqrt(np.mean(dG**2)), max(np.abs(dG)), np.sqrt(np.mean(dGfrac**2)), max(np.abs(dGfrac)))
    return G

#    shutil.rmtree(odir)
            if not os.path.exists(odir): os.makedirs(odir)
            apath = os.path.join(odir, "current")
            submit_psi(apath, d, mvals)
            for p in range(self.FF.np):
                def subjob(mvals_,h):
                    apath = os.path.join(odir, str(p), str(h))
                    submit_psi(apath, d, mvals_)
                    #logger.info("Will set up a job for %s, parameter %i\n" % (d, p))
                    return 0.0
                if self.callderivs[d][p]:
                    if AHess:
                        f12d3p(fdwrap(subjob, mvals, p, h=self.h), h = self.h, f0 = 0.0)
                    elif AGrad:
                        if self.bidirect:
                            f12d3p(fdwrap(subjob, mvals, p, h=self.h), h = self.h, f0 = 0.0)
                        else:
                            f1d2p(fdwrap(subjob, mvals, p, h=self.h), h = self.h, f0 = 0.0)

GDx      = np.zeros((FF.np,length))
    GDy      = np.zeros((FF.np,length))
    GDz      = np.zeros((FF.np,length))
    if not AGrad:
        return G, GDx, GDy, GDz
    def energy_driver(mvals_):
        FF.make(mvals_)
        if dipole:
            return engine.energy_dipole()
        else:
            return engine.energy()

    ED0      = energy_driver(mvals)
    for i in pgrad:
        logger.info("%i %s\r" % (i, (FF.plist[i] + " "*30)))
        EDG, _   = f12d3p(fdwrap(energy_driver,mvals,i),h,f0=ED0)
        if dipole:
            G[i,:]   = EDG[:,0]
            GDx[i,:] = EDG[:,1]
            GDy[i,:] = EDG[:,2]
            GDz[i,:] = EDG[:,3]
        else:
            G[i,:]   = EDG[:]
    return G, GDx, GDy, GDz

QBN = np.dot(self.qmboltz_wts[:NS],self.boltz_wts[:NS])
        #==============================================================#
        #             STEP 2: Loop through the snapshots.              #
        #==============================================================#
        if self.all_at_once:
            logger.debug("\rExecuting\r")
            M_all = self.energy_force_transform()
            if self.asym:
                M_all[:, 0] -= M_all[self.smin, 0]
            if not cv and (AGrad or AHess):
                def callM(mvals_):
                    logger.debug("\r")
                    pvals = self.FF.make(mvals_)
                    return self.energy_force_transform()
                for p in self.pgrad:
                    dM_all[:,p,:], ddM_all[:,p,:] = f12d3p(fdwrap(callM, mvals, p), h = self.h, f0 = M_all)
                    if self.asym:
                        dM_all[:, p, 0] -= dM_all[self.smin, p, 0]
                        ddM_all[:, p, 0] -= ddM_all[self.smin, p, 0]
        if self.force and not in_fd():
            self.maxfatom = -1
            self.maxfshot = -1
            self.maxdf = 0.0
        for i in range(NS):
            if i % 100 == 0:
                logger.debug("\rIncrementing quantities for snapshot %i\r" % i)
            # Build Boltzmann weights and increment partition function.
            P   = self.boltz_wts[i]
            Z  += P
            R   = self.qmboltz_wts[i]*self.boltz_wts[i] / QBN
            Y  += R
            # Recall reference (QM) data

"""
    def single_point(mvals_):
        FF.make(mvals_)
        if dipole:
            return engine.energy_dipole()
        else:
            return engine.energy()

    ED0 = single_point(mvals)
    G   = OrderedDict()
    G['potential'] = np.zeros((FF.np, ED0.shape[0]))
    if dipole:
        G['dipole'] = np.zeros((FF.np, ED0.shape[0], 3))
    for i in pgrad:
        logger.info("%i %s\r" % (i, (FF.plist[i] + " "*30)))
        edg, _ = f12d3p(fdwrap(single_point,mvals,i),h,f0=ED0)
        if dipole:
            G['potential'][i] = edg[:,0]
            G['dipole'][i]    = edg[:,1:]
        else:
            G['potential'][i] = edg[:]
    return G