Python forcefield.FF類代碼示例

def main():
    ## Set some basic options.  Note that 'forcefield' requires 'ffdir'
    ## which indicates the relative path of the force field.
    options, tgt_opts = parse_inputs(argv[1])
    MyFF = FF(options)
    Prec=int(argv[2])
    if 'read_mvals' in options:
        mvals = np.array(options['read_mvals'])
    else:
        mvals = np.zeros(len(MyFF.pvals0))
    MyFF.make(mvals,False,'NewFF',precision=Prec)

def energy_derivatives(engine, FF, mvals, h, pgrad, length, AGrad=True, dipole=False):

    """
    Compute the first and second derivatives of a set of snapshot
    energies with respect to the force field parameters.

    This basically calls the finite difference subroutine on the
    energy_driver subroutine also in this script.

    @param[in] mvals Mathematical parameter values
    @param[in] h Finite difference step size
    @param[in] phase The phase (liquid, gas) to perform the calculation on
    @param[in] AGrad Switch to turn derivatives on or off; if off, return all zeros
    @param[in] dipole Switch for dipole derivatives.
    @return G First derivative of the energies in a N_param x N_coord array
    @return GDx First derivative of the box dipole moment x-component in a N_param x N_coord array
    @return GDy First derivative of the box dipole moment y-component in a N_param x N_coord array
    @return GDz First derivative of the box dipole moment z-component in a N_param x N_coord array

    """
    G        = np.zeros((FF.np,length))
    GDx      = np.zeros((FF.np,length))
    GDy      = np.zeros((FF.np,length))
    GDz      = np.zeros((FF.np,length))
    if not AGrad:
        return G, GDx, GDy, GDz
    def energy_driver(mvals_):
        FF.make(mvals_)
        if dipole:
            return engine.energy_dipole()
        else:
            return engine.energy()

    ED0      = energy_driver(mvals)
    for i in pgrad:
        logger.info("%i %s\r" % (i, (FF.plist[i] + " "*30)))
        EDG, _   = f12d3p(fdwrap(energy_driver,mvals,i),h,f0=ED0)
        if dipole:
            G[i,:]   = EDG[:,0]
            GDx[i,:] = EDG[:,1]
            GDy[i,:] = EDG[:,2]
            GDz[i,:] = EDG[:,3]
        else:
            G[i,:]   = EDG[:]
    #reset FF parameters
    FF.make(mvals)
    return G, GDx, GDy, GDz

    def runTest(self):
        """Check implicit hydration free energy study (Hydration target) converges to expected results"""
        self.logger.debug("\nSetting input file to 'optimize.in'\n")
        input_file='optimize.in'

        ## The general options and target options that come from parsing the input file
        self.logger.debug("Parsing inputs...\n")
        options, tgt_opts = parse_inputs(input_file)
        self.logger.debug("options:\n%s\n\ntgt_opts:\n%s\n\n" % (str(options), str(tgt_opts)))

        self.assertEqual(dict,type(options), msg="\nParser gave incorrect type for options")
        self.assertEqual(list,type(tgt_opts), msg="\nParser gave incorrect type for tgt_opts")
        for target in tgt_opts:
            self.assertEqual(dict, type(target), msg="\nParser gave incorrect type for target dict")

        ## The force field component of the project
        self.logger.debug("Creating forcefield using loaded options: ")
        forcefield  = FF(options)
        self.logger.debug(str(forcefield) + "\n")
        self.assertEqual(FF, type(forcefield), msg="\nExpected forcebalance forcefield object")

        ## The objective function
        self.logger.debug("Creating object using loaded options and forcefield: ")
        objective   = Objective(options, tgt_opts, forcefield)
        self.logger.debug(str(objective) + "\n")
        self.assertEqual(Objective, type(objective), msg="\nExpected forcebalance objective object")

        ## The optimizer component of the project
        self.logger.debug("Creating optimizer: ")
        optimizer   = Optimizer(options, objective, forcefield)
        self.logger.debug(str(optimizer) + "\n")
        self.assertEqual(Optimizer, type(optimizer), msg="\nExpected forcebalance optimizer object")

        ## Actually run the optimizer.
        self.logger.debug("Done setting up! Running optimizer...\n")
        result = optimizer.Run()

        self.logger.debug("\nOptimizer finished. Final results:\n")
        self.logger.debug(str(result) + '\n')

        self.assertNdArrayEqual(EXPECTED_ETHANOL_RESULTS,forcefield.create_pvals(result),delta=0.02,
                                msg="\nCalculation results have changed from previously calculated values.\n"
                                "If this seems reasonable, update EXPECTED_ETHANOL_RESULTS in test_system.py with these values")

def energy_driver(mvals,FF,xyz,tky,verbose=False,dipole=False):
    """
    Compute a set of snapshot energies (and optionally, dipoles) as a function of the force field parameters.

    ForceBalance creates the force field, TINKER reads it in, and we loop through the snapshots
    to compute the energies.

    @param[in] mvals Mathematical parameter values
    @param[in] FF ForceBalance force field object
    @return E A numpy array of energies in kilojoules per mole

    """
    # Part of the command line argument to TINKER.
    basename = xyz[:-4]
    xin = "%s" % xyz + ("" if tky == None else " -k %s" % tky)
    xain = "%s.arc" % basename + ("" if tky == None else " -k %s" % tky)
    
    # Print the force field file from the ForceBalance object, with modified parameters.
    FF.make(mvals)
    
    # Execute TINKER.
    cmdstr = "./analyze %s" % xain
    oanl = _exec(cmdstr,stdin="E",print_command=verbose,print_to_screen=verbose)

    # Read potential energy from file.
    E = []
    for line in oanl:
        if 'Total Potential Energy : ' in line:
            E.append(float(line.split()[4]))
    E = np.array(E) * 4.184
    if dipole:
        # If desired, read dipole from file.
        D = []
        for line in oanl:
            if 'Dipole X,Y,Z-Components :' in line:
                D.append([float(line.split()[i]) for i in range(-3,0)])
        D = np.array(D)
        # Return a Nx4 array with energies in the first column and dipole in columns 2-4.
        answer = np.hstack((E.reshape(-1,1), D.reshape(-1,3)))
        return answer
    else:
        return E

def main():

    """
    Usage: (runcuda.sh) npt.py <openmm|gromacs|tinker> <liquid_nsteps> <liquid_timestep (fs)> <liquid_intvl (ps> <temperature> <pressure>

    This program is meant to be called automatically by ForceBalance on
    a GPU cluster (specifically, subroutines in openmmio.py).  It is
    not easy to use manually.  This is because the force field is read
    in from a ForceBalance 'FF' class.

    I wrote this program because automatic fitting of the density (or
    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.

    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(open('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)
    liquid_timestep = TgtOptions['liquid_timestep']
    liquid_nsteps = TgtOptions['liquid_md_steps']
    liquid_nequil = TgtOptions['liquid_eq_steps']
    liquid_intvl = TgtOptions['liquid_interval']
    liquid_fnm = TgtOptions['liquid_coords']
    gas_timestep = TgtOptions['gas_timestep']
    gas_nsteps = TgtOptions['gas_md_steps']
    gas_nequil = TgtOptions['gas_eq_steps']
    gas_intvl = TgtOptions['gas_interval']
    gas_fnm = TgtOptions['gas_coords']

    # Number of threads, multiple timestep integrator, anisotropic box etc.
    threads = TgtOptions.get('md_threads', 1)
    mts = TgtOptions.get('mts_integrator', 0)
    rpmd_beads = TgtOptions.get('rpmd_beads', 0)
    force_cuda = TgtOptions.get('force_cuda', 0)
    anisotropic = TgtOptions.get('anisotropic_box', 0)
    minimize = TgtOptions.get('minimize_energy', 1)

    # Print all options.
    printcool_dictionary(TgtOptions, title="Options from ForceBalance")
    liquid_snapshots = (liquid_nsteps * liquid_timestep / 1000) / liquid_intvl
    liquid_iframes = 1000 * liquid_intvl / liquid_timestep
    gas_snapshots = (gas_nsteps * gas_timestep / 1000) / gas_intvl
    gas_iframes = 1000 * gas_intvl / gas_timestep
    logger.info("For the condensed phase system, I will collect %i snapshots spaced apart by %i x %.3f fs time steps\n" \
        % (liquid_snapshots, liquid_iframes, liquid_timestep))
    if liquid_snapshots < 2:
        raise Exception('Please set the number of liquid time steps so that you collect at least two snapshots (minimum %i)' \
                            % (2000 * (liquid_intvl/liquid_timestep)))
    logger.info("For the gas phase system, I will collect %i snapshots spaced apart by %i x %.3f fs time steps\n" \
        % (gas_snapshots, gas_iframes, gas_timestep))
    if gas_snapshots < 2:
        raise Exception('Please set the number of gas time steps so that you collect at least two snapshots (minimum %i)' \
                            % (2000 * (gas_intvl/gas_timestep)))

    #----
    # Loading coordinates
    #----
    ML = Molecule(liquid_fnm)
    MG = Molecule(gas_fnm)
    # Determine the number of molecules in the condensed phase coordinate file.
    NMol = len(ML.molecules)

    #----
    # Setting up MD simulations
    #----
    EngOpts = OrderedDict()
    EngOpts["liquid"] = OrderedDict([("coords", liquid_fnm), ("mol", ML), ("pbc", True)])
    EngOpts["gas"] = OrderedDict([("coords", gas_fnm), ("mol", MG), ("pbc", False)])
    GenOpts = OrderedDict([('FF', FF)])
    if engname == "openmm":
        # OpenMM-specific options
        EngOpts["liquid"]["platname"] = 'CUDA'
        EngOpts["gas"]["platname"] = 'Reference'
#.........這裏部分代碼省略.........

 def energy_driver(mvals_):
     FF.make(mvals_)
     return engine.energy_dipole()

 def rpmd_energy_driver(mvals_):
     FF.make(mvals_)
     if dipole:
         return engine.energy_dipole_rpmd()
     else:
         return engine.energy_rpmd()

def main():

    """
    Usage: (runcuda.sh) nvt.py <openmm|gromacs|tinker> <liquid_nsteps> <liquid_timestep (fs)> <liquid_intvl (ps> <temperature>

    This program is meant to be called automatically by ForceBalance on
    a GPU cluster (specifically, subroutines in openmmio.py).  It is
    not easy to use manually.  This is because the force field is read
    in from a ForceBalance 'FF' class.
    """

    printcool("ForceBalance condensed phase NVT 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)
    nvt_timestep = TgtOptions['nvt_timestep']
    nvt_md_steps = TgtOptions['nvt_md_steps']
    nvt_eq_steps = TgtOptions['nvt_eq_steps']
    nvt_interval = TgtOptions['nvt_interval']
    liquid_fnm = TgtOptions['nvt_coords']

    # Number of threads, multiple timestep integrator, anisotropic box etc.
    threads = TgtOptions.get('md_threads', 1)
    mts = TgtOptions.get('mts_integrator', 0)
    rpmd_beads = TgtOptions.get('rpmd_beads', 0)
    force_cuda = TgtOptions.get('force_cuda', 0)
    nbarostat = TgtOptions.get('n_mcbarostat', 25)
    anisotropic = TgtOptions.get('anisotropic_box', 0)
    minimize = TgtOptions.get('minimize_energy', 1)

    # Print all options.
    printcool_dictionary(TgtOptions, title="Options from ForceBalance")
    nvt_snapshots = (nvt_timestep * nvt_md_steps / 1000) / nvt_interval
    nvt_iframes = 1000 * nvt_interval / nvt_timestep
    logger.info("For the condensed phase system, I will collect %i snapshots spaced apart by %i x %.3f fs time steps\n" \
        % (nvt_snapshots, nvt_iframes, nvt_timestep))
    if nvt_snapshots < 2:
        raise Exception('Please set the number of liquid time steps so that you collect at least two snapshots (minimum %i)' \
                            % (2000 * (nvt_interval/nvt_timestep)))

    #----
    # Loading coordinates
    #----
    ML = Molecule(liquid_fnm, toppbc=True)
    # Determine the number of molecules in the condensed phase coordinate file.
    NMol = len(ML.molecules)
    TgtOptions['n_molecules'] = NMol
    logger.info("There are %i molecules in the liquid\n" % (NMol))

    #----
    # Setting up MD simulations
    #----
    EngOpts = OrderedDict()
    EngOpts["liquid"] = OrderedDict([("coords", liquid_fnm), ("mol", ML), ("pbc", True)])
    if "nonbonded_cutoff" in TgtOptions:
        EngOpts["liquid"]["nonbonded_cutoff"] = TgtOptions["nonbonded_cutoff"]
    else:
        largest_available_cutoff = min(ML.boxes[0][:3]) / 2 - 0.1
        EngOpts["liquid"]["nonbonded_cutoff"] = largest_available_cutoff
        logger.info("nonbonded_cutoff is by default set to the largest available value: %g Angstrom" %largest_available_cutoff)
    if "vdw_cutoff" in TgtOptions:
        EngOpts["liquid"]["vdw_cutoff"] = TgtOptions["vdw_cutoff"]
    # Hard Code nonbonded_cutoff to 13A for test
    #EngOpts["liquid"]["nonbonded_cutoff"] = EngOpts["liquid"]["vdw_cutoff"] = 13.0
    GenOpts = OrderedDict([('FF', FF)])
    if engname == "openmm":
        # OpenMM-specific options
        EngOpts["liquid"]["platname"] = TgtOptions.get("platname", 'CUDA')
        if force_cuda:
            try: Platform.getPlatformByName('CUDA')
            except: raise RuntimeError('Forcing failure because CUDA platform unavailable')
            EngOpts["liquid"]["platname"] = 'CUDA'
        if threads > 1: logger.warn("Setting the number of threads will have no effect on OpenMM engine.\n")

    EngOpts["liquid"].update(GenOpts)
    for i in EngOpts:
        printcool_dictionary(EngOpts[i], "Engine options for %s" % i)

    # Set up MD options
    MDOpts = OrderedDict()
    MDOpts["liquid"] = OrderedDict([("nsteps", nvt_md_steps), ("timestep", nvt_timestep),
                                    ("temperature", temperature),
                                    ("nequil", nvt_eq_steps), ("minimize", minimize),
#.........這裏部分代碼省略.........

def main():

    """
    Run the script with -h for help
    Usage: python npt_tinker.py input.xyz [-k input.key] liquid_production_steps liquid_timestep liquid_interval temperature(K) pressure(atm)
    """

    if not os.path.exists(args.liquid_xyzfile):
        warn_press_key("Warning: %s does not exist, script cannot continue" % args.liquid_xyzfile)

    # Set up some conversion factors
    # All units are in kJ/mol
    N = niterations
    # Conversion factor for kT derived from:
    # In [6]: 1.0 / ((1.0 * kelvin * BOLTZMANN_CONSTANT_kB * AVOGADRO_CONSTANT_NA) / kilojoule_per_mole)
    # Out[6]: 120.27221251395186
    T     = temperature
    mBeta = -120.27221251395186 / temperature
    Beta  =  120.27221251395186 / temperature
    kT    =  0.0083144724712202 * temperature
    # Conversion factor for pV derived from:
    # In [14]: 1.0 * atmosphere * nanometer ** 3 * AVOGADRO_CONSTANT_NA / kilojoule_per_mole
    # Out[14]: 0.061019351687175
    pcon  =  0.061019351687175

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

    #=================================================================#
    #     Get the number of molecules from the liquid xyz file.       #
    #=================================================================#

    xin = "%s" % args.liquid_xyzfile + ("" if args.liquid_keyfile == None else " -k %s" % args.liquid_keyfile)
    cmdstr = "./analyze %s" % xin
    oanl = _exec(cmdstr,stdin="G",print_command=True,print_to_screen=True)
    molflag = False
    for line in oanl:
        if 'Number of Molecules' in line:
            if not molflag:
                NMol = int(line.split()[-1])
                molflag = True
            else:
                raise Exception("TINKER output contained more than one line with the words 'Number of Molecules'")
    if molflag:
        print "Detected %i Molecules" % NMol
    if not molflag:
        raise Exception("Failed to detect the number of molecules")

    #=================================================================#
    # Run the simulation for the full system and analyze the results. #
    #=================================================================#
    Rhos, Potentials, Kinetics, Volumes, Dips = run_simulation(args.liquid_xyzfile,args.liquid_keyfile,tstep=timestep,nstep=nsteps,neq=nequiliterations,npr=niterations,verbose=True)
    Energies = Potentials + Kinetics
    V  = Volumes
    pV = pressure * Volumes
    H = Energies + pV

    # Get the energy and dipole gradients.
    print "Post-processing the liquid simulation snapshots."
    G, GDx, GDy, GDz = energy_dipole_derivatives(mvals,h,FF,args.liquid_xyzfile,args.liquid_keyfile,AGrad)
    print

    #==============================================#
    # Now run the simulation for just the monomer. #
    #==============================================#
    _a, mPotentials, mKinetics, _b, _c = run_simulation(args.gas_xyzfile,args.gas_keyfile,tstep=m_timestep,nstep=m_nsteps,neq=m_nequiliterations,npr=m_niterations,pbc=False)
    mEnergies = mPotentials + mKinetics
    mN = len(mEnergies)
    print "Post-processing the gas simulation snapshots."
    mG = energy_derivatives(mvals,h,FF,args.gas_xyzfile,args.gas_keyfile,AGrad)
    print

    numboots = 1000    
    def bootstats(func,inputs):
        # Calculate error using bootstats method
        dboot = []
        for i in range(numboots):
            newins = {k : v[np.random.randint(len(v),size=len(v))] for k,v in inputs.items()}
            dboot.append(np.mean(func(**newins)))
        return func(**inputs),np.std(np.array(dboot))
        
    def calc_arr(b = None, **kwargs):
        # This tomfoolery is required because of Python syntax;
        # default arguments must come after nondefault arguments
        # and kwargs must come at the end.  This function is used
        # in bootstrap error calcs and also in derivative calcs.
        if 'arr' in kwargs:
            arr = kwargs['arr']
        if b == None: b = np.ones(len(arr),dtype=float)
        return bzavg(arr,b)

    # The density in kg/m^3.
    # Note: Not really necessary to use bootstrap here, but good to 
    # demonstrate the principle.
    Rho_avg,  Rho_err  = bootstats(calc_arr,{'arr':Rhos})
    Rho_err *= np.sqrt(statisticalInefficiency(Rhos))

#.........這裏部分代碼省略.........