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

print "Minimizing the energy... (starting energy % .3f kJ/mol)" % simulation.context.getState(getEnergy=True).getPotentialEnergy().value_in_unit(kilojoule_per_mole),
    simulation.minimizeEnergy()
    print "Done (final energy % .3f kJ/mol)" % simulation.context.getState(getEnergy=True).getPotentialEnergy().value_in_unit(kilojoule_per_mole)
    # Assign velocities.
    velocities = generateMaxwellBoltzmannVelocities(system, temperature)
    simulation.context.setVelocities(velocities)
    if verbose:
        # Print out the platform used by the context
        print "I'm using the platform", simulation.context.getPlatform().getName()
        # Print out the properties of the platform
        printcool_dictionary({i:simulation.context.getPlatform().getPropertyValue(simulation.context,i) for i in simulation.context.getPlatform().getPropertyNames()},title="Platform %s has properties:" % simulation.context.getPlatform().getName())
    # Serialize the system if we want.
    Serialize = 0
    if Serialize:
        serial = XmlSerializer.serializeSystem(system)
        with wopen('system.xml') as f: f.write(serial)
    #==========================================#
    # Computing a bunch of initial values here #
    #==========================================#
    if pbc:
        # Show initial system volume.
        box_vectors = system.getDefaultPeriodicBoxVectors()
        volume = compute_volume(box_vectors)
        if verbose: print "initial system volume = %.1f nm^3" % (volume / nanometers**3)
    # Determine number of degrees of freedom.
    kB = BOLTZMANN_CONSTANT_kB * AVOGADRO_CONSTANT_NA
    # The center of mass motion remover is also a constraint.
    ndof = 3*system.getNumParticles() - system.getNumConstraints() - 3
    # Compute total mass.
    mass = compute_mass(system).in_units_of(gram / mole) /  AVOGADRO_CONSTANT_NA # total system mass in g
    # Initialize statistics.
    data = dict()

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, pdb, FF, Xyzs, Settings, Sim, kT, calc_eps0, {'d_':Dips,'v_':V}, Boxes)
        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)

    with wopen(os.path.join('npt_result.p')) as f: lp_dump((Rhos, Volumes, Energies, Dips, G, [GDx, GDy, GDz], mEnergies, mG, Rho_err, Hvap_err, Alpha_err, Kappa_err, Cp_err, Eps0_err, NMol),f)

#throw error if more than one script
        #write script into .txt file and parse as text
        fflist = list(self.ffdata[ffname].iter())
        scriptElements = [elem for elem in fflist if elem.tag=='Script']
        if len(scriptElements) > 1:
            logger.error('XML file'+ffname+'contains more than one script! Consolidate your scripts into one script!\n')
            raise RuntimeError
        elif len(scriptElements)==1:
            Script = scriptElements[0].text
            ffnameList = ffname.split('.')
            ffnameScript = ffnameList[0]+'Script.txt'
            absScript = os.path.join(self.root, self.ffdir, ffnameScript)
            if os.path.exists(absScript):
                logger.error('XML file '+absScript+' already exists on disk! Please delete it\n')
                raise RuntimeError
            wfile = forcebalance.nifty.wopen(absScript)
            wfile.write(Script)
            wfile.close()
            self.addff(ffnameScript, xmlScript=True)
            os.unlink(absScript)

        for e in self.ffdata[ffname].getroot().xpath('//@parameterize/..'):
            parameters_to_optimize = [i.strip() for i in e.get('parameterize').split(',')]
            for p in parameters_to_optimize:
                if p not in e.attrib:
                    logger.error("Parameter \'%s\' is not found for \'%s\', please check %s" % (p, e.get('type'), ffname) )
                    raise RuntimeError
                pid = self.Readers[ffname].build_pid(e, p)
                self.map[pid] = self.np
                # offxml file later than v0.3 may have unit strings in the field
                quantity_str = e.get(p)
                res = re.search(r'^[-+]?[0-9]*\.?[0-9]*([eEdD][-+]?[0-9]+)?', quantity_str)

line_out.append(line)
    # Sanity checks: If frcmod and mol2 files are provided to this function,
    # they should be in the leap.cmd file as well.  There should be exactly
    # one PDB file being loaded.
    for i in frcmod:
        if i not in have_fmod:
            warn_press_key("WARNING: %s is not being loaded in %s" % (i, fnm))
    for i in mol2:
        if i not in have_mol2:
            warn_press_key("WARNING: %s is not being loaded in %s" % (i, fnm))

    fout = fnm+'_'
    line_out.append('saveamberparm %s %s.prmtop %s.inpcrd\n' % (ambername, prefix, prefix))
    line_out.append('quit\n')
    if os.path.exists(fout): os.remove(fout)
    with wopen(fout) as f: print(''.join(line_out), file=f)

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)

    with wopen(os.path.join('npt_result.p')) as f: 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),f)

# Average over all points
        Objective += np.sum(Objs, axis=0)
        Gradient  += np.sum(Grads, axis=0)
        Hessian   += np.sum(Hess, axis=0)
        
        # Store gradients and setup print map 
        GradMapPrint = [["#Point"] + self.FF.plist]

        for pt in self.points:
            temp  = pt.temperature
            press = pt.pressure
            GradMapPrint.append([' %8.2f %8.1f' % (temp, press)] +
                                ["% 9.3e" % i for i in grads[pt.idnr-1]])

        o = wopen('gradient_%s.dat' % quantity)
        for line in GradMapPrint:
            print >> o, ' '.join(line)
        o.close()
        
        printer = OrderedDict([("    %-5d %-12.2f %-8.1f"
                                % (pt.idnr, pt.temperature, pt.pressure),
                                ("% -10.3f % -10.3f  +- %-8.3f % -8.3f % -9.5f % -9.5f"
                                 % (Exp[pt.idnr-1], values[pt.idnr-1],
                                    errors[pt.idnr-1], Deltas[pt.idnr-1],
                                    Weights[pt.idnr-1], Objs[pt.idnr-1])))
                                    for pt in self.points])
                
        return { "X": Objective, "G": Gradient, "H": Hessian, "info": printer }

# Create instances of the MD Engine objects.
    Engine = EngineClass(name=Sim.type, **Sim.EngOpts)
    click() # Start timer.
    # This line runs the condensed phase simulation.
    Engine.molecular_dynamics(**Sim.MDOpts)
    logger.info("MD simulation took %.3f seconds\n" % click())
    # Extract properties.
    Results = Engine.md_extract(OrderedDict([(i, {}) for i in Sim.timeseries.keys()]))
    # Calculate energy and dipole derivatives if needed.
    if AGrad:
        Results['derivatives'] = energy_derivatives(Engine, FF, mvals, h, pgrad, dipole='dipole' in Sim.timeseries.keys())
    # Dump results to file
    logger.info("Writing final force field.\n")
    pvals = FF.make(mvals)
    logger.info("Writing all simulation data to disk.\n")
    with wopen('md_result.p') as f:
        lp_dump(Results, f)

def writechk(self):
        """ Write the checkpoint file for the main optimizer. """
        if self.wchk_fnm is not None:
            logger.info("Writing the checkpoint file %s\n" % self.wchk_fnm)
            with wopen(os.path.join(self.root,self.wchk_fnm), binary=True) as f: pickle.dump(self.chk, f)