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

print "Creating Double Precision Simulation for parameter derivatives"
        mSim, _ = create_simulation_object(mpdb, mSettings, pbc=False, precision="double")
    mG = energy_derivatives(mvals, h, mpdb, FF, mXyzs, mSettings, mSim, None, AGrad)

    # pV_avg and mean(pV) are exactly the same.
    pV = (pressure * Data['volume'] * AVOGADRO_CONSTANT_NA).value_in_unit(kilojoule_per_mole)
    kT = (kB * temperature).value_in_unit(kilojoule_per_mole)

    # The enthalpy of vaporization in kJ/mol.
    Hvap_avg = mPot_avg - Pot_avg / 216 + kT - np.mean(pV) / 216
    Hvap_err = np.sqrt(Pot_err**2 / 216**2 + mPot_err**2 + pV_err**2/216**2)

    # 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)) / 216

    print "The finite difference step size is:",h

    Sep = printcool("Density: % .4f +- % .4f kg/m^3, Analytic Derivative" % (Rho_avg, Rho_err))
    FF.print_map(vals=GRho)
    print Sep

    H = Energies + pV
    V = np.array(Volumes)
    numboots = 1000
    L = len(H)

# Run the simulation for the full system and analyze the results. #
   #=================================================================#
   Data, Xyzs, Boxes, Rhos, Energies = run_simulation(pdb, direct_kwargs)
   # Get statistics from our simulation.
   Rho_avg, Rho_err, Pot_avg, Pot_err = analyze(Data)
   # Now that we have the coordinates, we can compute the energy derivatives.
   G, Hd = energy_derivatives(mvals, h, pdb, FF, Xyzs, direct_kwargs, Boxes)
   # The density derivative can be computed using the energy derivative.
   N = len(Xyzs)
   kB = units.BOLTZMANN_CONSTANT_kB * units.AVOGADRO_CONSTANT_NA
   mBeta = (-1 / (temperature * kB)).value_in_unit(units.mole / units.kilojoule)
   Beta = (1 / (temperature * kB)).value_in_unit(units.mole / units.kilojoule)
   # It becomes necessary to make the gradient squared for the second derivative.
   GG = G * G
   # Build the first density derivative .
   GRho = mBeta * (flat(np.mat(G) * col(Rhos)) / N - np.mean(Rhos) * np.mean(G, axis=1))
   # The second density derivative is so inaccurate we might as well not compute it. -_-
   # HdRho = mBeta * (flat(np.mat(Hd) * col(Rhos)) / N
   #                  - Beta * flat(np.mat(GG) * col(Rhos)) / N
   #                  + Beta * (flat(np.mat(G) * col(Rhos)) / N) * np.mean(G, axis=1)
   #                  - np.mean(Rhos) * (np.mean(Hd, axis=1)
   #                                     - Beta * np.mean(GG, axis=1)
   #                                     + Beta * np.mean(G, axis=1) * np.mean(G, axis=1)))

   #==============================================#
   # Now run the simulation for just the monomer. #
   #==============================================#
   mpdb = PDBFile('mono.pdb')
   mData, mXyzs, _trash, _crap, mEnergies = run_simulation(mpdb, mono_kwargs, pbc=False, Trajectory=False)
   # Get statistics from our simulation.
   _trash, _crap, mPot_avg, mPot_err = analyze(mData)
   # Now that we have the coordinates, we can compute the energy derivatives.

Rho_avg, Rho_err, Pot_avg, Pot_err, pV_avg, pV_err = analyze(Data)
    # Now that we have the coordinates, we can compute the energy derivatives.
    # First create a double-precision simulation object.
    DoublePrecisionDerivatives = True
    if DoublePrecisionDerivatives and AGrad:
        print "Creating Double Precision Simulation for parameter derivatives"
        Sim, _ = create_simulation_object(pdb, Settings, pbc=True, precision="double")
    G, GDx, GDy, GDz = energy_dipole_derivatives(mvals, h, pdb, FF, Xyzs, Settings, Sim, Boxes, AGrad)
    # The density derivative can be computed using the energy derivative.
    N = len(Xyzs)
    kB = BOLTZMANN_CONSTANT_kB * AVOGADRO_CONSTANT_NA
    T = temperature / kelvin
    mBeta = (-1 / (temperature * kB)).value_in_unit(mole / kilojoule)
    Beta = (1 / (temperature * kB)).value_in_unit(mole / kilojoule)
    # Build the first density derivative .
    GRho = mBeta * (flat(np.mat(G) * col(Rhos)) / N - np.mean(Rhos) * np.mean(G, axis=1))

    #==============================================#
    # Now run the simulation for just the monomer. #
    #==============================================#
    global timestep, nsteps, niterations, nequiliterations
    timestep = 1.0 * femtosecond       # timestep for integration
    nsteps   = 100                     # number of steps per data record
    nequiliterations = 500             # number of equilibration iterations
    niterations = 5000                # number of iterations to collect data for

    mpdb = PDBFile('mono.pdb')
    mData, mXyzs, _trash, _crap, mEnergies, _nah, _dontneed, mSim = run_simulation(mpdb, mSettings, pbc=False, Trajectory=False)
    # Get statistics from our simulation.
    _trash, _crap, mPot_avg, mPot_err, __trash, __crap = analyze(mData)
    # Now that we have the coordinates, we can compute the energy derivatives.
    if DoublePrecisionDerivatives and AGrad:

h_ = kwargs['h_']
        if 'v_' in kwargs:
            v_ = kwargs['v_']
        return 1/(kT*T) * (bzavg(h_*v_,b)-bzavg(h_,b)*bzavg(v_,b))/bzavg(v_,b)
    Alpha = calc_alpha(None, **{'h_':H, 'v_':V})
    Alphaboot = []
    for i in range(numboots):
        boot = np.random.randint(L,size=L)
        Alphaboot.append(calc_alpha(None, **{'h_':H[boot], 'v_':V[boot]}))
    Alphaboot = np.array(Alphaboot)
    Alpha_err = np.std(Alphaboot) * max([np.sqrt(statisticalInefficiency(V)),np.sqrt(statisticalInefficiency(H))])

    ## Thermal expansion coefficient analytic derivative
    GAlpha1 = mBeta * covde(H*V) / avg(V)
    GAlpha2 = Beta * avg(H*V) * covde(V) / avg(V)**2
    GAlpha3 = flat(np.mat(G)*col(V))/N/avg(V) - Gbar
    GAlpha4 = Beta * covde(H)
    GAlpha  = (GAlpha1 + GAlpha2 + GAlpha3 + GAlpha4)/(kT*T)
    Sep = printcool("Thermal expansion coefficient: % .4e +- %.4e K^-1\nAnalytic Derivative:" % (Alpha, Alpha_err))
    FF.print_map(vals=GAlpha)
    if FDCheck:
        GAlpha_fd = property_derivatives(mvals, h, pdb, FF, Xyzs, Settings, Sim, kT, calc_alpha, {'h_':H,'v_':V}, Boxes)
        Sep = printcool("Numerical Derivative:")
        FF.print_map(vals=GAlpha_fd)
        Sep = printcool("Difference (Absolute, Fractional):")
        absfrac = ["% .4e  % .4e" % (i-j, (i-j)/j) for i,j in zip(GAlpha, GAlpha_fd)]
        FF.print_map(vals=absfrac)

    ## Isothermal compressibility
    bar_unit = 1.0*bar*nanometer**3/kilojoules_per_mole/item
    def calc_kappa(b=None, **kwargs):
        if b == None: b = np.ones(L,dtype=float)

dVdqP   = np.matrix(self.invdists[i])
            espqval = espqvals[i]
            espmval = dVdqP * col(new_charges(mvals))
            desp    = flat(espmval) - espqval
            X      += P * np.dot(desp, desp) / self.nesp
            Q      += P * np.dot(espqval, espqval) / self.nesp
            D      += P * (np.dot(espqval, espqval) / self.nesp - (np.sum(espqval) / self.nesp)**2)
            if AGrad:
                dVdqM   = (dVdqP * dqPdqM).T
                for p, vsd in ddVdqPdVS.items():
                    dVdqM[p,:] += flat(vsd[i] * col(new_charges(mvals)))
                G      += flat(P * 2 * dVdqM * col(desp)) / self.nesp
                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

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)))
            
        Answer = {'X':X,'G':G,'H':H}
        return Answer

if self.backup:
                            for fnm in self.FF.fnms:
                                if os.path.exists(os.path.join(self.resdir, fnm)):
                                    bak(os.path.join(self.resdir, fnm))
                        self.FF.make(xk,printdir=self.resdir)
                        # logger.info("The force field has been written to the '%s' directory.\n" % self.resdir)
                        outfnm = self.save_mvals_to_input(xk)
                        # logger.info("Input file with optimization parameters saved to %s.\n" % outfnm)
                        printcool("The force field has been written to the %s directory.\n"
                                  "Input file with optimization parameters saved to %s." % (self.resdir, outfnm), color=0)
                    #================================#
                    #|  Hessian update for BFGS.    |#
                    #================================#
                    if b_BFGS:
                        Hnew = H_stor.copy()
                        Dx   = col(xk - xk_prev)
                        Dy   = col(G  - G_prev)
                        Mat1 = (np.dot(Dy,Dy.T))/(np.dot(Dy.T,Dx))[0,0]
                        Mat2 = (np.dot(np.dot(Hnew,Dx),np.dot(Hnew,Dx).T))/(multi_dot([Dx.T,Hnew,Dx]))[0,0]
                        Hnew += Mat1-Mat2
                        H = Hnew.copy()
                        data['H'] = H.copy()
                    # (Experimental): Deleting lines in the parameter file
                    if len(self.FF.prmdestroy_this) > 0:
                        self.FF.prmdestroy_save.append(self.FF.prmdestroy_this)
                        self.FF.linedestroy_save.append(self.FF.linedestroy_this) 
            # Update objective function history.
            X_hist = np.append(X_hist, X)
            # Take the stdev over the previous (hist) values.
            # Multiply by 2, so when hist=2 this is simply the difference.
            stdfront = np.std(X_hist[-self.hist:]) if len(X_hist) > self.hist else np.std(X_hist)
            stdfront *= 2

dz = d_[:,2]
        D2  = bzavg(dx**2,b)-bzavg(dx,b)**2
        D2 += bzavg(dy**2,b)-bzavg(dy,b)**2
        D2 += bzavg(dz**2,b)-bzavg(dz,b)**2
        return prefactor*D2/bzavg(v_,b)/T

    Eps0, Eps0_err = bootstats(calc_eps0,{'d_':Dips, 'v_':V})
    Eps0 += 1.0
    Eps0_err *= 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.mat(GDx)*col(Dx))/N - avg(Dx)*(np.mean(GDx,axis=1))) - Beta*(covde(Dx**2) - 2*avg(Dx)*covde(Dx))
    GD2 += 2*(flat(np.mat(GDy)*col(Dy))/N - avg(Dy)*(np.mean(GDy,axis=1))) - Beta*(covde(Dy**2) - 2*avg(Dy)*covde(Dy))
    GD2 += 2*(flat(np.mat(GDz)*col(Dz))/N - 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(mvals, h, FF, args.liquid_xyzfile, args.liquid_keyfile, 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)

    ## Print the final force field.
    pvals = FF.make(mvals)