How to use the forcebalance.nifty.flat function in forcebalance

else:
            # G0 and H0 are used for determining the expected function change.
            G0 = G.copy()
            H0 = H.copy()
            G = np.delete(G, self.excision)
            H = np.delete(H, self.excision, axis=0)
            H = np.delete(H, self.excision, axis=1)
            
            logger.debug("Inverting Hessian:\n")                 ###
            logger.debug(" G:\n")                                ###
            pvec1d(G,precision=5, loglevel=DEBUG)                ###
            logger.debug(" H:\n")                                ###
            pmat2d(H,precision=5, loglevel=DEBUG)                ###
            
            Hi = invert_svd(np.mat(H))
            dx = flat(-1 * Hi * col(G))
            
            logger.debug(" dx:\n")                               ###
            pvec1d(dx,precision=5, loglevel=DEBUG)                     ###
            # dxa = -solve(H, G)          # Take Newton Raphson Step ; use -1*G if want steepest descent.
            # dxa = flat(dxa)
            # print " dxa:"                              ###
            # pvec1d(dxa,precision=5)                    ###
            
            logger.info('\n')                                      ###
            for i in self.excision:    # Reinsert deleted coordinates - don't take a step in those directions
                dx = np.insert(dx, i, 0)
            def para_solver(L):
                # Levenberg-Marquardt
                # HT = H + (L-1)**2*np.diag(np.diag(H))
                # Attempt to use plain Levenberg
                HT = H + (L-1)**2*np.eye(len(H))

def para_solver(L):
                # Levenberg-Marquardt
                # HT = H + (L-1)**2*np.diag(np.diag(H))
                # Attempt to use plain Levenberg
                HT = H + (L-1)**2*np.eye(len(H))
                logger.debug("Inverting Scaled Hessian:\n")
                logger.debug(" G:\n")
                pvec1d(G,precision=5, loglevel=DEBUG)
                logger.debug(" HT: (Scal = %.4f)\n" % (1+(L-1)**2))
                pmat2d(HT,precision=5, loglevel=DEBUG)
                Hi = invert_svd(HT)
                dx = flat(-1 * np.dot(Hi, col(G)))
                logger.debug(" dx:\n")
                pvec1d(dx,precision=5, loglevel=DEBUG)
                sol = flat(0.5*multi_dot([row(dx), H, col(dx)]))[0] + np.dot(dx,G)
                for i in self.excision:    # Reinsert deleted coordinates - don't take a step in those directions
                    dx = np.insert(dx, i, 0)
                return dx, sol

# Build the first Hvap derivative.
   # We don't pass it back, but nice for printing.
   GHvap = np.mean(G,axis=1)
   GHvap += mBeta * (flat(np.mat(G) * col(Energies)) / N - Pot_avg * np.mean(G, axis=1))
   GHvap /= 216
   GHvap -= np.mean(mG,axis=1)
   GHvap -= mBeta * (flat(np.mat(mG) * col(mEnergies)) / N - mPot_avg * np.mean(mG, axis=1))
   GHvap *= -1
   
   GHvap -= mBeta * (flat(np.mat(G) * col(pV)) / N - np.mean(pV) * np.mean(G, axis=1))

   print "pV_avg = ", pV_avg, "pV_err = ", pV_err
   print "kT, -pV", kT, -1 * np.mean(pV)
   print "Derivative of pV term (note it is SUBTRACTED from enthalpy of vaporization, so the contribution to Hvap is minus one times this):"
   FF.print_map(vals = mBeta * (flat(np.mat(G) * col(pV)) / N - np.mean(pV) * np.mean(G, axis=1)))

   bar = printcool("Density: % .4f +- % .4f kg/m^3, Derivatives below" % (Rho_avg, Rho_err))
   FF.print_map(vals=GRho)
   print bar
   
   print "Box energy:", np.mean(Energies)
   print "Monomer energy:", np.mean(mEnergies)
   bar = printcool("Enthalpy of Vaporization: % .4f +- %.4f kJ/mol, Derivatives below" % (Hvap_avg, Hvap_err))
   FF.print_map(vals=GHvap)
   print bar
   # Print the final force field.
   pvals = FF.make(mvals)

   with open(os.path.join('npt_result.p'),'w') as f: lp_dump((Rhos, pV, Energies, G, mEnergies, mG, Rho_err, Hvap_err),f)

Eps0 = calc_eps0(None,**{'d_':Dips, 'v_':V})
    Eps0boot = []
    for i in range(numboots):
        boot = np.random.randint(L,size=L)
        Eps0boot.append(calc_eps0(None,**{'d_':Dips[boot], 'v_':V[boot]}))
    Eps0boot = np.array(Eps0boot)
    Eps0_err = np.std(Eps0boot)*np.sqrt(np.mean([statisticalInefficiency(Dips[:,0]),statisticalInefficiency(Dips[:,1]),statisticalInefficiency(Dips[:,2])]))
 
    # Dielectric constant analytic derivative
    Dx = Dips[:,0]
    Dy = Dips[:,1]
    Dz = Dips[:,2]
    D2 = avg(Dx**2)+avg(Dy**2)+avg(Dz**2)-avg(Dx)**2-avg(Dy)**2-avg(Dz)**2
    GD2  = 2*(flat(np.dot(GDx,col(Dx)))/L - avg(Dx)*(np.mean(GDx,axis=1))) - Beta*(covde(Dx**2) - 2*avg(Dx)*covde(Dx))
    GD2 += 2*(flat(np.dot(GDy,col(Dy)))/L - avg(Dy)*(np.mean(GDy,axis=1))) - Beta*(covde(Dy**2) - 2*avg(Dy)*covde(Dy))
    GD2 += 2*(flat(np.dot(GDz,col(Dz)))/L - avg(Dz)*(np.mean(GDz,axis=1))) - Beta*(covde(Dz**2) - 2*avg(Dz)*covde(Dz))
    GEps0 = prefactor*(GD2/avg(V) - mBeta*covde(V)*D2/avg(V)**2)/T
    Sep = printcool("Dielectric constant:           % .4e +- %.4e\nAnalytic Derivative:" % (Eps0, Eps0_err))
    FF.print_map(vals=GEps0)
    if FDCheck:
        GEps0_fd = property_derivatives(Liquid, FF, mvals, h, pgrad, kT, calc_eps0, {'d_':Dips,'v_':V})
        Sep = printcool("Numerical Derivative:")
        FF.print_map(vals=GEps0_fd)
        Sep = printcool("Difference (Absolute, Fractional):")
        absfrac = ["% .4e  % .4e" % (i-j, (i-j)/j) for i,j in zip(GEps0,GEps0_fd)]
        FF.print_map(vals=absfrac)

    logger.info("Writing final force field.\n")
    pvals = FF.make(mvals)

    logger.info("Writing all simulation data to disk.\n")
    lp_dump((Rhos, Volumes, Potentials, Energies, Dips, G, [GDx, GDy, GDz], mPotentials, mEnergies, mG, Rho_err, Hvap_err, Alpha_err, Kappa_err, Cp_err, Eps0_err, NMol),'npt_result.p')

dEg = results['Potential_Derivatives']
                        dEaq = results['Potential_Derivatives'] + results['Hydration_Derivatives']
                        data[p]['dHyd'] = (flat(np.dot(dEaq,col(expmbH))/L)-np.mean(dEg,axis=1)*np.mean(expmbH)) / np.mean(expmbH)
                elif p == "liq":
                    Eg = results['Potentials'] - results['Hydration']
                    Eaq = results['Potentials']
                    # Mean and standard error of the exponentiated hydration energy.
                    exppbH = np.exp(+1.0*beta*results['Hydration'])
                    data[p]['Hyd'] = +kT*np.log(np.mean(exppbH))
                    # Estimate standard error by bootstrap method.  We also multiply by the 
                    # square root of the statistical inefficiency of the hydration energy time series.
                    data[p]['HydErr'] = np.std([+kT*np.log(np.mean(exppbH[np.random.randint(L,size=L)])) for i in range(100)]) * np.sqrt(statisticalInefficiency(results['Hydration']))
                    if AGrad: 
                        dEg = results['Potential_Derivatives'] - results['Hydration_Derivatives']
                        dEaq = results['Potential_Derivatives']
                        data[p]['dHyd'] = -(flat(np.dot(dEg, col(exppbH))/L)-np.mean(dEaq,axis=1)*np.mean(exppbH)) / np.mean(exppbH)
                os.chdir('..')
            # Calculate the hydration free energy using gas phase, liquid phase or the average of both.
            # Note that the molecular dynamics methods return energies in kJ/mol.
            if self.hfemode == 'exp_gas':
                self.hfe_dict[label] = data['gas']['Hyd'] / 4.184
                self.hfe_err[label] = data['gas']['HydErr'] / 4.184
            elif self.hfemode == 'exp_liq':
                self.hfe_dict[label] = data['liq']['Hyd'] / 4.184
                self.hfe_err[label] = data['liq']['HydErr'] / 4.184
            elif self.hfemode == 'exp_both':
                self.hfe_dict[label] = 0.5*(data['liq']['Hyd']+data['gas']['Hyd']) / 4.184
                self.hfe_err[label] = 0.5*(data['liq']['HydErr']+data['gas']['HydErr']) / 4.184
            if AGrad:
                # Calculate the derivative of the hydration free energy.
                if self.hfemode == 'exp_gas':
                    dD[:, ilabel] = self.whfe[ilabel]*data['gas']['dHyd'] / 4.184

G = zeros(np, dtype=float)
        H = zeros((np, np), dtype=float)
        for i in range(self.ns):
            P   = self.whamboltz_wts[i]
            Z  += P
            dVdqP   = mat(self.invdists[i])
            espqval = espqvals[i]
            espmval = dVdqP * col(getqatoms(mvals))
            desp    = flat(espmval) - espqval
            X      += P * dot(desp, desp) / self.nesp
            D      += P * (dot(espqval, espqval) / self.nesp - (sum(espqval) / self.nesp)**2)
            if AGrad:
                dVdqM   = (dVdqP * dqPdqM).T
                for p, vsd in ddVdqPdVS.items():
                    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

dddVdqPdVS2 = {}
        if AGrad:
            for p in range(np):
                if 'VSITE' in self.FF.plist[p]:
                    ddVdqPdVS[p], dddVdqPdVS2[p] = f12d3p(fdwrap(self.build_invdist,mvals,p), h = self.h, f0 = self.invdists)
        X = 0
        D = 0
        G = zeros(np, dtype=float)
        H = zeros((np, np), dtype=float)
        for i in range(self.ns):
            P   = self.whamboltz_wts[i]
            Z  += P
            dVdqP   = mat(self.invdists[i])
            espqval = espqvals[i]
            espmval = dVdqP * col(getqatoms(mvals))
            desp    = flat(espmval) - espqval
            X      += P * dot(desp, desp) / self.nesp
            D      += P * (dot(espqval, espqval) / self.nesp - (sum(espqval) / self.nesp)**2)
            if AGrad:
                dVdqM   = (dVdqP * dqPdqM).T
                for p, vsd in ddVdqPdVS.items():
                    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

def para_solver(L):
                # Levenberg-Marquardt
                # HT = H + (L-1)**2*np.diag(np.diag(H))
                # Attempt to use plain Levenberg
                HT = H + (L-1)**2*np.eye(len(H))
                logger.debug("Inverting Scaled Hessian:\n")
                logger.debug(" G:\n")
                pvec1d(G,precision=5, loglevel=DEBUG)
                logger.debug(" HT: (Scal = %.4f)\n" % (1+(L-1)**2))
                pmat2d(HT,precision=5, loglevel=DEBUG)
                Hi = invert_svd(HT)
                dx = flat(-1 * np.dot(Hi, col(G)))
                logger.debug(" dx:\n")
                pvec1d(dx,precision=5, loglevel=DEBUG)
                sol = flat(0.5*multi_dot([row(dx), H, col(dx)]))[0] + np.dot(dx,G)
                for i in self.excision:    # Reinsert deleted coordinates - don't take a step in those directions
                    dx = np.insert(dx, i, 0)
                return dx, sol

def _compute(self, dx):
                    self.dx = dx.copy()
                    Tmp = np.mat(self.H)*col(dx)
                    Reg_Term   = self.Penalty.compute(xkd+flat(dx), Obj0)
                    self.Val   = (X + np.dot(dx, G) + 0.5*row(dx)*Tmp + Reg_Term[0] - data['X'])[0,0]
                    self.Grad  = flat(col(G) + Tmp) + Reg_Term[1]
                def compute_val(self, dx):

HT = H + (L-1)**2*np.eye(len(H))
                logger.debug("Inverting Scaled Hessian:\n")                       ###
                logger.debug(" G:\n")                                             ###
                pvec1d(G,precision=5, loglevel=DEBUG)                                   ###
                logger.debug(" HT: (Scal = %.4f)\n" % (1+(L-1)**2))               ###
                pmat2d(HT,precision=5, loglevel=DEBUG)                                  ###
                Hi = invert_svd(np.mat(HT))
                dx = flat(-1 * Hi * col(G))
                logger.debug(" dx:\n")                                            ###
                pvec1d(dx,precision=5, loglevel=DEBUG)                                  ###
                # dxa = -solve(HT, G)
                # dxa = flat(dxa)
                # print " dxa:"                                           ###
                # pvec1d(dxa,precision=5)                                 ###
                # print                                                   ###
                sol = flat(0.5*row(dx)*np.mat(H)*col(dx))[0] + np.dot(dx,G)
                for i in self.excision:    # Reinsert deleted coordinates - don't take a step in those directions
                    dx = np.insert(dx, i, 0)
                return dx, sol