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

Of course this script can also be executed locally.  Just make sure
   you have the pickled 'forcebalance.p' file.

   """
   
   # Select platform.
   platform = openmm.Platform.getPlatformByName('Cuda')
   # Set the CUDA device to the environment variable or zero otherwise
   cuda_device = os.environ.get('CUDA_DEVICE',"0")
   print "Setting CUDA Device to", cuda_device
   platform.setPropertyDefaultValue("CudaDevice", cuda_device)
   # Specify the PDB here so we may make the Simulation class.
   pdb = PDBFile(sys.argv[1])
   # Load the force field in from the ForceBalance pickle.
   FF,mvals,h = lp_load(open('forcebalance.p'))
   # Create the force field XML files.
   pvals = FF.make(os.getcwd(),mvals,False)
   # Run the simulation.
   Data, Xyzs, Boxes, Rhos = run_simulation(pdb)

   #print Data['potential'].value_in_unit(units.kilojoule / units.mole)
   #energy_driver(mvals, pdb, FF, Xyzs, Boxes, verbose=True)
   # Get statistics from our simulation.
   analyze(Data)
   # Now that we have the coordinates, we can compute the energy derivatives.
   G, Hd = energy_derivatives(mvals, h, pdb, FF, Xyzs, 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)

calculations need to be distributed to the queue.  The main instance
   of ForceBalance (running on my workstation) queues up a bunch of these
   jobs (using Work Queue).  Then, I submit a bunch of workers to GPU
   clusters (e.g. Certainty, Keeneland).  The worker scripts connect to
   the main instance and receives one of these jobs.

   This script can also be executed locally, if you want to (e.g. for
   debugging).  Just make sure you have the pickled 'forcebalance.p'
   file.

   """
   
   # Create an OpenMM PDB object so we may make the Simulation class.
   pdb = PDBFile(sys.argv[1])
   # Load the force field in from the ForceBalance pickle.
   FF,mvals,h,AGrad = lp_load(open('forcebalance.p'))
   # Create the force field XML files.
   FF.make(mvals)

   #=================================================================#
   # Run the simulation for the full system and analyze the results. #
   #=================================================================#
   Data, Xyzs, Boxes, Rhos, Energies = run_simulation(pdb, mutual_kwargs if FF.amoeba_pol == 'mutual' else direct_kwargs, Trajectory=True)
   # Get statistics from our simulation.
   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.
   G, Hd = energy_derivatives(mvals, h, pdb, FF, Xyzs, mutual_kwargs if FF.amoeba_pol == 'mutual' else direct_kwargs, Boxes, AGrad)
   # 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)

other equilibrium properties) is computationally intensive, and the
   calculations need to be distributed to the queue.  The main instance
   of ForceBalance (running on my workstation) queues up a bunch of these
   jobs (using Work Queue).  Then, I submit a bunch of workers to GPU
   clusters (e.g. Certainty, Keeneland).  The worker scripts connect to
   the main instance and receives one of these jobs.

   Of course this script can also be executed locally.  Just make sure
   you have the pickled 'forcebalance.p' file.

   """
   
   # Specify the PDB here so we may make the Simulation class.
   pdb = PDBFile(sys.argv[1])
   # Load the force field in from the ForceBalance pickle.
   FF,mvals,h = lp_load(open('forcebalance.p'))
   # Create the force field XML files.
   pvals = FF.make(mvals,False)

   #=================================================================#
   # 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)

def get_exp(self, mvals, AGrad=False, AHess=False):
        """ Get the hydration free energy using the Zwanzig formula.  We will obtain two different estimates along with their uncertainties. """
        self.hfe_dict = OrderedDict()
        self.hfe_err = OrderedDict()
        dD = np.zeros((self.FF.np,len(self.IDs)))
        kT = (kb * self.hfe_temperature)
        beta = 1. / (kb * self.hfe_temperature)
        for ilabel, label in enumerate(self.IDs):
            os.chdir(label)
            # This dictionary contains observables keyed by each phase.
            data = defaultdict(dict)
            for p in ['gas', 'liq']:
                os.chdir(p)
                # Load the results from molecular dynamics.
                results = lp_load('md_result.p')
                L = len(results['Potentials'])
                if p == "gas":
                    Eg = results['Potentials']
                    Eaq = results['Potentials'] + results['Hydration']
                    # Mean and standard error of the exponentiated hydration energy.
                    expmbH = np.exp(-1.0*beta*results['Hydration'])
                    data[p]['Hyd'] = -kT*np.log(np.mean(expmbH))
                    # 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(expmbH[np.random.randint(L,size=L)])) for i in range(100)]) * np.sqrt(statisticalInefficiency(results['Hydration']))
                    if AGrad: 
                        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']

Parameters
        ----------
        dp : Point
            Store the calculated quantities in this point.

        Returns
        -------
        Nothing
        
        """
        abspath = os.path.join(os.getcwd(), '%d/md_result.p' % dp.idnr)

        if os.path.exists(abspath):
            logger.info('Reading data from ' + abspath + '.\n')

            vals, errs, grads = lp_load(abspath)

            dp.data["values"] = vals
            dp.data["errors"] = errs
            dp.data["grads"]  = grads

        else:
            msg = 'The file ' + abspath + ' does not exist so we cannot read it.\n'
            logger.warning(msg)

            dp.data["values"] = np.zeros((len(self.quantities)))
            dp.data["errors"] = np.zeros((len(self.quantities)))
            dp.data["grads"]  = np.zeros((len(self.quantities), self.FF.np))

from __future__ import print_function
import forcebalance
import forcebalance.objective
import forcebalance.nifty
from forcebalance.nifty import wopen
import tarfile
import os
import forcebalance.output
logger = forcebalance.output.getLogger("forcebalance")
logger.setLevel(forcebalance.output.DEBUG)

# load pickled variables from forcebalance.p
if os.path.exists('forcebalance.p'):
    mvals, AGrad, AHess, id_string, options, tgt_opts, forcefield, pgrad = forcebalance.nifty.lp_load('forcebalance.p')
else:
    forcefield, mvals = forcebalance.nifty.lp_load('forcefield.p')
    AGrad, AHess, id_string, options, tgt_opts, pgrad = forcebalance.nifty.lp_load('options.p')

print("Evaluating remote target ID: %s" % id_string)

options['root'] = os.getcwd()
forcefield.root = os.getcwd()

# set up forcefield
forcefield.make(mvals, printdir="forcefield")

# set up and evaluate target
tar = tarfile.open("target.tar.bz2", "r")
tar.extractall()
tar.close()

This script can also be executed locally, if you want to (e.g. for
    debugging).  Just make sure you have the pickled 'forcebalance.p'
    file.

    """

    printcool("ForceBalance condensed phase simulation using engine: %s" % engname.upper(), color=4, bold=True)

    #----
    # Load the ForceBalance pickle file which contains:
    #----
    # - Force field object
    # - Optimization parameters
    # - Options from the Target object that launched this simulation
    # - Switch for whether to evaluate analytic derivatives.
    FF,mvals,TgtOptions,AGrad = lp_load('forcebalance.p')
    FF.ffdir = '.'
    # Write the force field file.
    FF.make(mvals)

    #----
    # Load the options that are set in the ForceBalance input file.
    #----
    # Finite difference step size
    h = TgtOptions['h']
    pgrad = TgtOptions['pgrad']
    # MD options; time step (fs), production steps, equilibration steps, interval for saving data (ps)
    lipid_timestep = TgtOptions['lipid_timestep']
    lipid_nsteps = TgtOptions['lipid_md_steps']
    lipid_nequil = TgtOptions['lipid_eq_steps']
    lipid_intvl = TgtOptions['lipid_interval']
    lipid_fnm = TgtOptions['lipid_coords']

#!/usr/bin/env python

from mslib import MSMS
from forcebalance.nifty import lp_load, lp_dump
import numpy as np
import os

# Designed to be called from GenerateQMData.py
# I wrote this because MSMS seems to have a memory leak

xyz, radii, density = lp_load(open('msms_input.p'))
MS = MSMS(coords = list(xyz), radii = radii)
MS.compute(density=density)
vfloat, vint, tri = MS.getTriangles()
with open(os.path.join('msms_output.p'), 'w') as f: lp_dump(vfloat, f)

parser.add_argument('-min', '--minimize', dest='minimize', action='store_true',
                    help='Whether to minimize the energy before starting the simulation')
parser.add_argument('-o', '-out', '--output', dest='output', type=str, nargs='+', 
                    help='Specify the time series which are written to disk')

# Parse the command line options and save as a dictionary (don't save NoneTypes)
parsed = parser.parse_args()
args = OrderedDict([(i, j) for i, j in vars(parsed).items() if j is not None])

#----
# Load the ForceBalance pickle file which contains:
#----
# - Force field object
# - Optimization parameters
# - Options loaded from file
FF, mvals = lp_load('forcefield.p')
#----
# Load the simulation pickle file which contains:
#----
# - Target options
# - Engine options
# - MD simulation options
TgtOpts, EngOpts, MDOpts = lp_load('simulation.p')
FF.ffdir = '.'

# Engine name.
engname = TgtOpts['engname']

# Import modules and create the correct Engine object.
if engname == "openmm":
    try:
        from simtk.unit import *

parser.add_argument('-min', '--minimize', dest='minimize', action='store_true',
                    help='Whether to minimize the energy before starting the simulation')
parser.add_argument('-o', '-out', '--output', dest='output', type=str, nargs='+', 
                    help='Specify the time series which are written to disk')

# Parse the command line options and save as a dictionary (don't save NoneTypes)
parsed = parser.parse_args()
args = OrderedDict([(i, j) for i, j in vars(parsed).items() if j != None])

#----
# Load the ForceBalance pickle file which contains:
#----
# - Force field object
# - Optimization parameters
# - Options loaded from file
FF, mvals, Sim = lp_load(open('forcebalance.p'))
FF.ffdir = '.'

# Engine name.
engname = Sim.engname

# Import modules and create the correct Engine object.
if engname == "openmm":
    try:
        from simtk.unit import *
        from simtk.openmm import *
        from simtk.openmm.app import *
    except:
        traceback.print_exc()
        raise Exception("Cannot import OpenMM modules")
    from forcebalance.openmmio import *
    EngineClass = OpenMM