How to use forcebalance - 10 common examples

for line in config:
            line = line.strip()
            if line:
                line = line.split('=')
                headless_options[line[0]]=line[1]

    os.chdir(os.path.dirname(__file__) + "/..")

    if not os.path.exists('/tmp/forcebalance'): os.mkdir('/tmp/forcebalance')
    warningHandler = forcebalance.output.CleanFileHandler('/tmp/forcebalance/test.err','w')
    warningHandler.setLevel(forcebalance.output.WARNING)
    logfile = "/tmp/forcebalance/%s.log" % time.strftime('%m%d%y_%H%M%S')
    debugHandler = forcebalance.output.CleanFileHandler(logfile,'w')
    debugHandler.setLevel(forcebalance.output.DEBUG)
    
    forcebalance.output.getLogger("forcebalance.test").addHandler(warningHandler)
    forcebalance.output.getLogger("forcebalance.test").addHandler(debugHandler)

    options['loglevel']=forcebalance.output.DEBUG
    
    runner=ForceBalanceTestRunner()
    results=runner.run(**options)

    if headless_options.has_key('enable_smtp')\
    and headless_options['enable_smtp'].lower() in ['true','error']:
        if headless_options['enable_smtp'].lower()=='true' or not results.wasSuccessful():
            import smtplib
            from email.mime.text import MIMEText
            from email.mime.multipart import MIMEMultipart

            # establish connection with smtp server
            server = smtplib.SMTP(host = headless_options["smtp_server"],

def run_simulation(pdb,settings,pbc=True,Trajectory=True):
    """ Run a NPT simulation and gather statistics. """
    simulation, system = create_simulation_object(pdb, settings, pbc, "single")
    # Set initial positions.
    simulation.context.setPositions(pdb.positions)
    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.

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)

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()

# for i in range(numboots):
    #    boot = np.random.randint(L,size=L)
    #    Rhoboot.append(calc_rho(None,**{'r_':Rhos[boot]}))
    # Rhoboot = np.array(Rhoboot)
    # Rho_err = np.std(Rhoboot)
    if FDCheck:
        Sep = printcool("Numerical Derivative:")
        GRho1 = property_derivatives(mvals, h, pdb, FF, Xyzs, Settings, Sim, kT, calc_rho, {'r_':Rhos}, Boxes)
        FF.print_map(vals=GRho1)
        Sep = printcool("Difference (Absolute, Fractional):")
        absfrac = ["% .4e  % .4e" % (i-j, (i-j)/j) for i,j in zip(GRho, GRho1)]
        FF.print_map(vals=absfrac)

    print "Box energy:", np.mean(Energies)
    print "Monomer energy:", np.mean(mEnergies)
    Sep = printcool("Enthalpy of Vaporization: % .4f +- %.4f kJ/mol, Derivatives below" % (Hvap_avg, Hvap_err))
    FF.print_map(vals=GHvap)
    print Sep

    # Define some things to make the analytic derivatives easier.
    Gbar = np.mean(G,axis=1)
    def deprod(vec):
        return flat(np.mat(G)*col(vec))/N
    def covde(vec):
        return flat(np.mat(G)*col(vec))/N - Gbar*np.mean(vec)
    def avg(vec):
        return np.mean(vec)

    ## Thermal expansion coefficient and bootstrap error estimation
    def calc_alpha(b = None, **kwargs):
        if b == None: b = np.ones(L,dtype=float)
        if 'h_' in kwargs:

def test_fd_stencils(self):
        """Check finite difference stencils return approximately correct results"""
        func = lambda x: x[0]**2
        fd_stencils = [function for function in dir(forcebalance.finite_difference) if re.match('^f..?d.p$',function)]

        self.logger.debug("Comparing fd stencils against some simple functions\n")
        for func in self.functions:
            for p in range(1):
                for x in range(10):
                    input = [0,0,0]
                    input[p]=x
                    f=forcebalance.finite_difference.fdwrap(func[0], input, p)
                    for stencil in fd_stencils:
                        fd = eval("forcebalance.finite_difference.%s" % stencil)
                        result = fd(f,.0001)
                        if re.match('^f..d.p$', stencil):
                            self.assertAlmostEqual(result[0], func[1](input,p), places=3)
                            self.assertAlmostEqual(result[1], func[2](input,p), places=3)
                        else:
                            self.assertAlmostEqual(result, func[1](input,p), places=3)

def __init__(self, logger=forcebalance.output.getLogger("forcebalance.test"), verbose = False):
        self.logger = logger
        forcebalance.output.getLogger("forcebalance.test")

def __init__(self):
        """Add logging capabilities to the standard TestResult implementation"""
        super(ForceBalanceTestResult,self).__init__()
        self.logger = forcebalance.output.getLogger('forcebalance.test.results')

self.logger.debug("Verify nifty matrix manipulations perform as expected\n")
        ##matrix manipulations
        X=flat(X)
        self.assertEqual(X.shape, (6,))
        X=row(X)
        self.assertEqual(X.shape, (1,6))
        X=col(X)
        self.assertEqual(X.shape, (6,1))

        self.logger.debug("Running some test processes using nifty._exec()\n")
        ##_exec
        self.assertEqual(type(_exec("")),list)
        self.assertEqual(_exec("echo test")[0],"test")
        _exec("touch .test")
        self.assertTrue(os.path.isfile(".test"))
        _exec("rm .test")
        self.assertFalse(os.path.isfile(".test"))
        self.assertRaises(Exception, _exec, "exit 255")

    @param[in] pdb OpenMM PDB object
    @param[in] FF ForceBalance force field object
    @param[in] xyzs List of OpenMM positions
    @param[in] settings OpenMM settings for creating the System
    @param[in] boxes Periodic box vectors
    @return G First derivative of the energies in a N_param x N_coord array

    """

    G        = np.zeros((FF.np,len(xyzs)))
    if not AGrad:
        return G
    E0       = energy_driver(mvals, pdb, FF, xyzs, settings, simulation, boxes)
    CheckFDPts = False
    for i in range(FF.np):
        G[i,:], _ = f12d3p(fdwrap(energy_driver,mvals,i,key=None,pdb=pdb,FF=FF,xyzs=xyzs,settings=settings,simulation=simulation,boxes=boxes),h,f0=E0)
        if CheckFDPts:
            # Check whether the number of finite difference points is sufficient.  Forward difference still gives residual error of a few percent.
            G1 = f1d7p(fdwrap(energy_driver,mvals,i,key=None,pdb=pdb,FF=FF,xyzs=xyzs,settings=settings,simulation=simulation,boxes=boxes),h)
            dG = G1 - G[i,:]
            dGfrac = (G1 - G[i,:]) / G[i,:]
            print "Parameter %3i 7-pt vs. central derivative : RMS, Max error (fractional) = % .4e % .4e (% .4e % .4e)" % (i, np.sqrt(np.mean(dG**2)), max(np.abs(dG)), np.sqrt(np.mean(dGfrac**2)), max(np.abs(dGfrac)))
    return G