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

# This ensures that the OrderedDict doesn't get out of order.
            for Tgt in self.Targets:
                self.ObjDict[Tgt.name] = None
            # Loop through the targets and compute the objective function for ones that are finished.
            while len(Need2Evaluate) > 0:
                for Tgt in Need2Evaluate:
                    if Tgt.wq_complete():
                        # List of functions that I can call.
                        Funcs   = [Tgt.get_X, Tgt.get_G, Tgt.get_H]
                        # Call the appropriate function
                        Ans = Funcs[Order](mvals, customdir=customdir)
                        # Print out the qualitative indicators
                        if verbose:
                            Tgt.meta_indicate(customdir=customdir)
                        # Note that no matter which order of function we call, we still increment the objective / gradient / Hessian the same way.
                        if not in_fd():
                            self.ObjDict[Tgt.name] = {'w' : Tgt.weight/self.WTot , 'x' : Ans['X']}
                        for i in range(3):
                            Objective[Letters[i]] += Ans[Letters[i]]*Tgt.weight/self.WTot
                        Need2Evaluate.remove(Tgt)
                        break
                    else:
                        pass
        else:
            wq = getWorkQueue()
            if wq is not None:
                wq_wait(wq)
            for Tgt in self.Targets:
                # The first call is always done at the midpoint.
                Tgt.bSave = True
                # List of functions that I can call.
                Funcs   = [Tgt.get_X, Tgt.get_G, Tgt.get_H]

else:
            wq = getWorkQueue()
            if wq is not None:
                wq_wait(wq)
            for Tgt in self.Targets:
                # The first call is always done at the midpoint.
                Tgt.bSave = True
                # List of functions that I can call.
                Funcs   = [Tgt.get_X, Tgt.get_G, Tgt.get_H]
                # Call the appropriate function
                Ans = Funcs[Order](mvals, customdir=customdir)
                # Print out the qualitative indicators
                if verbose:
                    Tgt.meta_indicate(customdir=customdir)
                # Note that no matter which order of function we call, we still increment the objective / gradient / Hessian the same way.
                if not in_fd():
                    self.ObjDict[Tgt.name] = {'w' : Tgt.weight/self.WTot , 'x' : Ans['X']}
                for i in range(3):
                    Objective[Letters[i]] += Ans[Letters[i]]*Tgt.weight/self.WTot
        # The target has evaluated at least once.
        for Tgt in self.Targets:
            Tgt.evaluated = True
        # Safeguard to make sure we don't have exact zeros on Hessian diagonal
        for i in range(self.FF.np):
            if Objective['H'][i,i] == 0.0:
                Objective['H'][i,i] = 1.0
        return Objective

if AHess:
                    d2VdqM2 = np.zeros(dVdqM.shape)
                    for p, vsd in dddVdqPdVS2.items():
                        d2VdqM2[p,:] += flat(vsd[i] * col(new_charges(mvals)))
                    H      += np.array(P * 2 * (dVdqM * dVdqM.T + d2VdqM2 * col(desp))) / self.nesp
        # Redundant but we keep it anyway
        D /= Z
        X /= Z
        X /= D
        G /= Z
        G /= D
        H /= Z
        H /= D
        Q /= Z
        Q /= D
        if not in_fd():
            self.esp_err = np.sqrt(X)
            self.esp_ref = np.sqrt(Q)
            self.esp_err_pct = self.esp_err / self.esp_ref
            
        # Following is the restraint part
        # RESP hyperbola "strength" parameter; 0.0005 is weak, 0.001 is strong
        # RESP hyperbola "tightness" parameter; don't need to change this
        a = self.resp_a
        b = self.resp_b
        q = new_charges(mvals)
        R   = a*np.sum((q**2 + b**2)**0.5 - b)
        dR  = a*q*(q**2 + b**2)**-0.5
        ddR = a*b**2*(q**2 + b**2)**-1.5
        self.respterm = R
        X += R
        if AGrad:

diff_bond = ((vref_bonds - vtar_bonds) * scale_bond).tolist() if n_bonds > 0 else []
                # objective contribution from angles
                vtar_angles = v_ic['angles']
                diff_angle = (periodic_diff(vref_angles, vtar_angles, 360) * scale_angle).tolist() if n_angles > 0 else []
                # objective contribution from dihedrals
                vtar_dihedrals = v_ic['dihedrals']
                diff_dihedral = (periodic_diff(vref_dihedrals, vtar_dihedrals, 360) * scale_dihedral).tolist() if n_dihedrals > 0 else []
                # objective contribution from improper dihedrals
                vtar_impropers = v_ic['impropers']
                diff_improper = (periodic_diff(vref_impropers, vtar_impropers, 360) * scale_improper).tolist() if n_impropers > 0 else []
                # combine objective values into a big result list
                sys_obj_list = diff_bond + diff_angle + diff_dihedral + diff_improper
                # extend the result v_obj_list by individual terms in this system
                v_obj_list += sys_obj_list
                # save print string
                if not in_fd():
                    # For printing, we group the RMSD by type
                    rmsd_bond = compute_rmsd(vref_bonds, vtar_bonds)
                    rmsd_angle = compute_rmsd(vref_angles, vtar_angles, v_periodic=360)
                    rmsd_dihedral = compute_rmsd(vref_dihedrals, vtar_dihedrals, v_periodic=360)
                    rmsd_improper = compute_rmsd(vref_impropers, vtar_impropers, v_periodic=360)
                    obj_total = sum(v**2 for v in sys_obj_list)
                    self.PrintDict[sysname] = "% 9.3f % 7.2f % 9.3f % 7.2f % 9.3f % 7.2f % 9.3f % 7.2f %17.3f" % (rmsd_bond, \
                        bond_denom, rmsd_angle, angle_denom, rmsd_dihedral, dihedral_denom, rmsd_improper, improper_denom, obj_total)
            return np.array(v_obj_list, dtype=float)

calc_momvals = self.unpack_moments(calc_moments)

        D = calc_momvals - ref_momvals
        dV = np.zeros((self.FF.np,len(calc_momvals)))

        if AGrad or AHess:
            for p in self.pgrad:
                dV[p,:], _ = f12d3p(fdwrap(get_momvals, mvals, p), h = self.h, f0 = calc_momvals)
                
        Answer['X'] = np.dot(D,D)
        for p in self.pgrad:
            Answer['G'][p] = 2*np.dot(D, dV[p,:])
            for q in self.pgrad:
                Answer['H'][p,q] = 2*np.dot(dV[p,:], dV[q,:])

        if not in_fd():
            self.calc_moments = calc_moments
            self.objective = Answer['X']

        return Answer

dNfrac = MBP * dN_M / qN_M + QBP * dN_Q / qN_Q
            dT = MBP * dT_M / Z + QBP * dT_Q / Y
            qT = MBP * qT_M / Z + QBP * qT_Q / Y
            dTfrac = MBP * dT_M / qT_M + QBP * dT_Q / qT_Q
            # if cv:
            #     dNfrac = MBP * sqrt(mean(array([SPiXi[i,i]/QQ_M[i,i] for i in range(1+3*self.fitatoms, 1+3*(self.fitatoms+self.nnf))]))) + \
            #         QBP * sqrt(mean(array([SRiXi[i,i]/QQ_Q[i,i] for i in range(1+3*self.fitatoms, 1+3*(self.fitatoms+self.nnf))])))
            #     dTfrac = MBP * sqrt(mean(array([SPiXi[i,i]/QQ_M[i,i] for i in range(1+3*(self.fitatoms+self.nnf), 1+3*(self.fitatoms+self.nnf+self.ntq))]))) + \
            #         QBP * sqrt(mean(array([SRiXi[i,i]/QQ_Q[i,i] for i in range(1+3*(self.fitatoms+self.nnf), 1+3*(self.fitatoms+self.nnf+self.ntq))])))
            # else:
            #     dNfrac = MBP * sqrt(mean(array([SPiXi[i]/QQ_M[i] for i in range(1+3*self.fitatoms, 1+3*(self.fitatoms+self.nnf)) if abs(QQ_M[i]) > 1e-3]))) + \
            #         QBP * sqrt(mean(array([SRiXi[i]/QQ_Q[i] for i in range(1+3*self.fitatoms, 1+3*(self.fitatoms+self.nnf)) if abs(QQ_Q[i]) > 1e-3])))
            #     dTfrac = MBP * sqrt(mean(array([SPiXi[i]/QQ_M[i] for i in range(1+3*(self.fitatoms+self.nnf), 1+3*(self.fitatoms+self.nnf+self.ntq)) if abs(QQ_M[i]) > 1e-3]))) + \
            #         QBP * sqrt(mean(array([SRiXi[i]/QQ_Q[i] for i in range(1+3*(self.fitatoms+self.nnf), 1+3*(self.fitatoms+self.nnf+self.ntq)) if abs(QQ_Q[i]) > 1e-3])))
        # Save values to qualitative indicator if not inside finite difference code.
        if not in_fd():
            self.e_ref = MBP * sqrt(QQ_M[0]/Z - Q0_M[0]**2/Z/Z) + QBP * sqrt((QQ_Q[0]/Y - Q0_Q[0]**2/Y/Y))
            self.e_err = dE
            self.e_err_pct = dEfrac
            if self.force:
                self.f_ref = qF
                self.f_err = dF
                self.f_err_pct = dFfrac
            if self.use_nft:
                self.nf_ref = qN
                self.nf_err = dN
                self.nf_err_pct = dNfrac
                self.tq_ref = qT
                self.tq_err = dT
                self.tq_err_pct = dTfrac
            pvals = self.FF.make(mvals) # Write a force field that isn't perturbed by finite differences.
        if cv:

dVdqM[p,:] += flat(vsd[i] * col(getqatoms(mvals)))
                G      += flat(P * 2 * dVdqM * col(desp)) / self.nesp
                if AHess:
                    d2VdqM2 = zeros(dVdqM.shape, dtype=float)
                    for p, vsd in dddVdqPdVS2.items():
                        d2VdqM2[p,:] += flat(vsd[i] * col(getqatoms(mvals)))
                    H      += array(P * 2 * (dVdqM * dVdqM.T + d2VdqM2 * col(desp))) / self.nesp
        # Redundant but we keep it anyway
        D /= Z
        X /= Z
        X /= D
        G /= Z
        G /= D
        H /= Z
        H /= D
        if not in_fd():
            self.esp_err = sqrt(X)
        # Following is the restraint part
        # RESP hyperbola "strength" parameter; 0.0005 is weak, 0.001 is strong
        # RESP hyperbola "tightness" parameter; don't need to change this
        a = self.resp_a
        b = self.resp_b
        q = getqatoms(mvals)
        R   = a*sum((q**2 + b**2)**0.5 - b)
        dR  = a*q*(q**2 + b**2)**-0.5
        ddR = a*b**2*(q**2 + b**2)**-1.5
        self.respterm = R
        X += R
        if AGrad:
            G += flat(dqPdqM.T * col(dR))
            if AHess:
                H += diag(flat(dqPdqM.T * col(ddR)))

else:
                            if wq != None:
                                # Since the calculations are submitted as 3-point finite difference, this part of the code
                                # actually only reads from half of the completed calculations.
                                grad[p] = f1d2p(fdwrap(reader, mvals, p, h=self.h), h = self.h, f0 = x)
                            else:
                                grad[p] = f1d2p(fdwrap(self.driver, mvals, p, d=d), h = self.h, f0 = x)
                            
            self.objd[d] = x
            self.gradd[d] = grad
            self.hdiagd[d] = hdiag
            X += x
            G += grad
            #H += np.diag(hdiag)
            H += hess
        if not in_fd():
            self.objective = X
            self.objvals = self.objd
        # print self.objd
        # print self.gradd
        # print self.hdiagd
                    
        if float('Inf') in pvals:
            return {'X' : 1e10, 'G' : G, 'H' : H}
        return {'X' : X, 'G' : G, 'H' : H}

plot_interaction_qm_vs_mm(qm_data, mm_data, title="Interaction Energy "+self.name)

        # Do the finite difference derivative.
        if AGrad or AHess:
            for p in self.pgrad:
                dV[p,:], _ = f12d3p(fdwrap(callM, mvals, p), h = self.h, f0 = emm)
            # Create the force field one last time.
            pvals  = self.FF.make(mvals)

        Answer['X'] = np.dot(self.prefactor*D/self.divisor,D/self.divisor)
        for p in self.pgrad:
            Answer['G'][p] = 2*np.dot(self.prefactor*D/self.divisor, dV[p,:]/self.divisor)
            for q in self.pgrad:
                Answer['H'][p,q] = 2*np.dot(self.prefactor*dV[p,:]/self.divisor, dV[q,:]/self.divisor)

        if not in_fd():
            self.emm = emm
            self.objective = Answer['X']

        ## QYD: try to clean up OpenMM engine.simulation objects to free up GPU memory
        try:
            if self.engine.name == 'openmm':
                if hasattr(self.engine, 'simulation'): del self.engine.simulation
                if hasattr(self.engine, 'A'): del self.engine.A
                if hasattr(self.engine, 'B'): del self.engine.B
        except:
            pass

        return Answer

def compute(mvals_):
            # This function has automatically assigned variable names from the interaction master file
            # Thus, all variable names in here are protected using an underscore.
            self.FF.make(mvals_)
            VectorD_ = []
            for sys_ in self.sys_opts:
                Energy_, RMSD_ = self.system_driver(sys_)
                # Energies are stored in a dictionary.
                EnergyDict[sys_] = Energy_
                RMSDNrm_ = RMSD_ / self.rmsd_denom
                w_ = self.sys_opts[sys_]['rmsd_weight'] if 'rmsd_weight' in self.sys_opts[sys_] else 1.0
                VectorD_.append(np.sqrt(w_)*RMSDNrm_)
                if not in_fd() and RMSD_ != 0.0:
                    self.RMSDDict[sys_] = "% 9.3f % 12.5f" % (RMSD_, w_*RMSDNrm_**2)
            VectorE_ = []
            for inter_ in self.inter_opts:
                def encloseInDictionary(matchobj):
                    return 'EnergyDict["' + matchobj.group(0)+'"]'
                # Here we need to evaluate a mathematical expression of the stored variables in EnergyDict.
                # We start by enclosing every variable in EnergyDict[""] and then calling eval on it.
                evalExpr = re.sub('[A-Za-z_][A-Za-z0-9_]*', encloseInDictionary, self.inter_opts[inter_]['equation'])
                Calculated_ = eval(evalExpr)
                Reference_ = self.inter_opts[inter_]['reference_physical']
                Delta_ = Calculated_ - Reference_
                Denom_ = self.energy_denom
                if self.cauchy:
                    Divisor_ = np.sqrt(Denom_**2 + Reference_**2)
                elif self.attenuate:
                    if Reference_ < Denom_: