当前位置: 首页>>代码示例>>Python>>正文


Python FF.print_map方法代码示例

本文整理汇总了Python中forcebalance.forcefield.FF.print_map方法的典型用法代码示例。如果您正苦于以下问题:Python FF.print_map方法的具体用法?Python FF.print_map怎么用?Python FF.print_map使用的例子?那么恭喜您, 这里精选的方法代码示例或许可以为您提供帮助。您也可以进一步了解该方法所在forcebalance.forcefield.FF的用法示例。


在下文中一共展示了FF.print_map方法的3个代码示例,这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞,您的评价将有助于系统推荐出更棒的Python代码示例。

示例1: main

# 需要导入模块: from forcebalance.forcefield import FF [as 别名]
# 或者: from forcebalance.forcefield.FF import print_map [as 别名]

#.........这里部分代码省略.........
    Potentials = prop_return['Potentials']

    #============================================#
    #  Compute the potential energy derivatives. #
    #============================================#
    if AGrad:
        logger.info("Calculating potential energy derivatives with finite difference step size: %f\n" % h)
        # Switch for whether to compute the derivatives two different ways for consistency.
        FDCheck = False
        printcool("Condensed phase energy and dipole derivatives\nInitializing array to length %i" % len(Potentials), color=4, bold=True)
        click()
        G, GDx, GDy, GDz = energy_derivatives(Liquid, FF, mvals, h, pgrad, len(Potentials), AGrad, dipole=False)
        logger.info("Condensed phase energy derivatives took %.3f seconds\n" % click())

    #==============================================#
    #  Condensed phase properties and derivatives. #
    #==============================================#

    # Physical constants
    kB = 0.008314472471220214
    T = temperature
    kT = kB * T # Unit: kJ/mol

    #--- Surface Tension ----
    logger.info("Start Computing surface tension.\n")
    perturb_proportion = 0.0005
    box_vectors = np.array(Liquid.xyz_omms[0][1]/nanometer) # obtain original box vectors from first frame
    delta_S = np.sqrt(np.sum(np.cross(box_vectors[0], box_vectors[1])**2)) * perturb_proportion * 2 # unit: nm^2. *2 for 2 surfaces
    # perturb xy area +
    click()
    scale_x = scale_y = np.sqrt(1 + perturb_proportion)
    scale_z = 1.0 / (1+perturb_proportion) # keep the box volumn while changing the area of xy plane
    Liquid.scale_box(scale_x, scale_y, scale_z)
    logger.info("scale_box+ took %.3f seconds\n" %click())
    # Obtain energies and gradients
    Potentials_plus = Liquid.energy()
    logger.info("Calculation of energies for perturbed box+ took %.3f seconds\n" %click())
    if AGrad:
        G_plus, _, _, _ = energy_derivatives(Liquid, FF, mvals, h, pgrad, len(Potentials), AGrad, dipole=False)
        logger.info("Calculation of energy gradients for perturbed box+ took %.3f seconds\n" %click())
    # perturb xy area - ( Note: also need to cancel the previous scaling)
    scale_x = scale_y = np.sqrt(1 - perturb_proportion) * (1.0/scale_x)
    scale_z = 1.0 / (1-perturb_proportion) * (1.0/scale_z)
    Liquid.scale_box(scale_x, scale_y, scale_z)
    logger.info("scale_box- took %.3f seconds\n" %click())
    # Obtain energies and gradients
    Potentials_minus = Liquid.energy()
    logger.info("Calculation of energies for perturbed box- took %.3f seconds\n" %click())
    if AGrad:
        G_minus, _, _, _ = energy_derivatives(Liquid, FF, mvals, h, pgrad, len(Potentials), AGrad, dipole=False)
        logger.info("Calculation of energy gradients for perturbed box- took %.3f seconds\n" %click())
    # Compute surface tension
    dE_plus = Potentials_plus - Potentials # Unit: kJ/mol
    dE_minus = Potentials_minus - Potentials # Unit: kJ/mol
    prefactor = -0.5 * kT / delta_S / 6.0221409e-1 # Unit mJ m^-2
    # Following equation: γ = -kT/(2ΔS) * [ ln<exp(-ΔE+/kT)> - ln<exp(-ΔE-/kT)> ]
    #plus_avg, plus_err = mean_stderr(np.exp(-dE_plus/kT))
    #minus_avg, minus_err = mean_stderr(np.exp(-dE_minus/kT))
    #surf_ten = -0.5 * kT / delta_S * ( np.log(plus_avg) - np.log(minus_avg) ) / 6.0221409e-1 # convert to mJ m^-2
    #surf_ten_err = 0.5 * kT / delta_S * ( np.log(plus_avg+plus_err) - np.log(plus_avg-plus_err) + np.log(minus_avg+minus_err) - np.log(minus_avg-minus_err) ) / 6.0221409e-1
    exp_dE_plus = np.exp(-dE_plus/kT)
    exp_dE_minus = np.exp(-dE_minus/kT)
    surf_ten = prefactor * ( np.log(np.mean(exp_dE_plus)) - np.log(np.mean(exp_dE_minus)) )
    # Use bootstrap method to estimate the error
    num_frames = len(exp_dE_plus)
    numboots = 1000
    surf_ten_boots = np.zeros(numboots)
    for i in xrange(numboots):
        boots_ordering = np.random.randint(num_frames, size=num_frames)
        boots_exp_dE_plus = np.take(exp_dE_plus, boots_ordering)
        boots_exp_dE_minus = np.take(exp_dE_minus, boots_ordering)
        surf_ten_boots[i] = prefactor * ( np.log(np.mean(boots_exp_dE_plus)) - np.log(np.mean(boots_exp_dE_minus)) )
    surf_ten_err = np.std(surf_ten_boots) * np.sqrt(np.mean([statisticalInefficiency(exp_dE_plus), statisticalInefficiency(exp_dE_minus)]))

    printcool("Surface Tension:       % .4f +- %.4f mJ m^-2" % (surf_ten, surf_ten_err))
    # Analytic Gradient of surface tension
    # Formula:      β = 1/kT
    #           ∂γ/∂α = -kT/(2ΔS) * { 1/<exp(-βΔE+)> * [<-β ∂E+/∂α exp(-βΔE+)> - <-β ∂E/∂α><exp(-βΔE+)>]
    #                                -1/<exp(-βΔE-)> * [<-β ∂E-/∂α exp(-βΔE-)> - <-β ∂E/∂α><exp(-βΔE-)>] }
    n_params = len(mvals)
    G_surf_ten = np.zeros(n_params)
    if AGrad:
        beta = 1.0 / kT
        plus_denom = np.mean(np.exp(-beta*dE_plus))
        minus_denom = np.mean(np.exp(-beta*dE_minus))
        for param_i in xrange(n_params):
            plus_left = np.mean(-beta * G_plus[param_i] * np.exp(-beta*dE_plus))
            plus_right = np.mean(-beta * G[param_i]) * plus_denom
            minus_left = np.mean(-beta * G_minus[param_i] * np.exp(-beta*dE_minus))
            minus_right = np.mean(-beta * G[param_i]) * minus_denom
            G_surf_ten[param_i] = prefactor * (1.0/plus_denom*(plus_left-plus_right) - 1.0/minus_denom*(minus_left-minus_right))
        printcool("Analytic Derivatives:")
        FF.print_map(vals=G_surf_ten)

    logger.info("Writing final force field.\n")
    pvals = FF.make(mvals)

    logger.info("Writing all results to disk.\n")
    result_dict = {'surf_ten': surf_ten, 'surf_ten_err': surf_ten_err, 'G_surf_ten': G_surf_ten}
    lp_dump(result_dict, 'nvt_result.p')
开发者ID:falagara,项目名称:forcebalance,代码行数:104,代码来源:nvt.py

示例2: main

# 需要导入模块: from forcebalance.forcefield import FF [as 别名]
# 或者: from forcebalance.forcefield.FF import print_map [as 别名]

#.........这里部分代码省略.........
    PrintEDA(mEDA, 1)

    #============================================#
    #  Compute the potential energy derivatives. #
    #============================================#
    logger.info("Calculating potential energy derivatives with finite difference step size: %f\n" % h)
    # Switch for whether to compute the derivatives two different ways for consistency.
    FDCheck = False

    # Create a double-precision simulation object if desired (seems unnecessary).
    DoublePrecisionDerivatives = False
    if engname == "openmm" and DoublePrecisionDerivatives and AGrad:
        logger.info("Creating Double Precision Simulation for parameter derivatives\n")
        Liquid = Engine(name="liquid", openmm_precision="double", **EngOpts["liquid"])
        Gas = Engine(name="gas", openmm_precision="double", **EngOpts["gas"])

    # Compute the energy and dipole derivatives.
    printcool("Condensed phase energy and dipole derivatives\nInitializing array to length %i" % len(Energies), color=4, bold=True)
    G, GDx, GDy, GDz = energy_derivatives(Liquid, FF, mvals, h, pgrad, len(Energies), AGrad, dipole=True)
    printcool("Gas phase energy derivatives", color=4, bold=True)
    mG, _, __, ___ = energy_derivatives(Gas, FF, mvals, h, pgrad, len(mEnergies), AGrad, dipole=False)

    #==============================================#
    #  Condensed phase properties and derivatives. #
    #==============================================#

    #----
    # Density
    #----
    # Build the first density derivative.
    GRho = mBeta * (flat(np.mat(G) * col(Rhos)) / L - np.mean(Rhos) * np.mean(G, axis=1))
    # Print out the density and its derivative.
    Sep = printcool("Density: % .4f +- % .4f kg/m^3\nAnalytic Derivative:" % (Rho_avg, Rho_err))
    FF.print_map(vals=GRho)
    logger.info(Sep)

    def calc_rho(b = None, **kwargs):
        if b == None: b = np.ones(L,dtype=float)
        if 'r_' in kwargs:
            r_ = kwargs['r_']
        return bzavg(r_,b)

    # No need to calculate error using bootstrap, but here it is anyway
    # Rhoboot = []
    # for i in range(numboots):
    #    boot = np.random.randint(N,size=N)
    #    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(Liquid, FF, mvals, h, pgrad, kT, calc_rho, {'r_':Rhos})
        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)

    #----
    # Enthalpy of vaporization
    #----
    H = Energies + pV
    V = np.array(Volumes)

    # Print out the liquid enthalpy.
    logger.info("Liquid enthalpy: % .4f kJ/mol, stdev % .4f ; (% .4f from energy, % .4f from pV)\n" % 
开发者ID:rbharath,项目名称:forcebalance,代码行数:70,代码来源:npt.py

示例3: main

# 需要导入模块: from forcebalance.forcefield import FF [as 别名]
# 或者: from forcebalance.forcefield.FF import print_map [as 别名]

#.........这里部分代码省略.........
    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))

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

    # The first density derivative
    GRho = mBeta * (flat(np.mat(G) * col(Rhos)) / N - np.mean(Rhos) * np.mean(G, axis=1))

    FDCheck = False

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

    if FDCheck:
        Sep = printcool("Numerical Derivative:")
        GRho1 = property_derivatives(mvals, h, FF, args.liquid_xyzfile, args.liquid_keyfile, kT, calc_arr, {'arr':Rhos})
        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)

    # The enthalpy of vaporization in kJ/mol.
    Ene_avg,  Ene_err  = bootstats(calc_arr,{'arr':Energies})
    mEne_avg, mEne_err = bootstats(calc_arr,{'arr':mEnergies})
    pV_avg,   pV_err   = bootstats(calc_arr,{'arr':pV})
    Ene_err  *= np.sqrt(statisticalInefficiency(Energies))
    mEne_err *= np.sqrt(statisticalInefficiency(mEnergies))
    pV_err   *= np.sqrt(statisticalInefficiency(pV))

    Hvap_avg = mEne_avg - Ene_avg / NMol + kT - np.mean(pV) / NMol
    Hvap_err = np.sqrt(Ene_err**2 / NMol**2 + mEne_err**2 + pV_err**2/NMol**2)

    # Build the first Hvap derivative.
    GHvap = np.mean(G,axis=1)
    GHvap += mBeta * (flat(np.mat(G) * col(Energies)) / N - Ene_avg * np.mean(G, axis=1))
    GHvap /= NMol
    GHvap -= np.mean(mG,axis=1)
    GHvap -= mBeta * (flat(np.mat(mG) * col(mEnergies)) / N - mEne_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)) / NMol

    print "Box total energy:", np.mean(Energies)
    print "Monomer total energy:", np.mean(mEnergies)
开发者ID:mlaury,项目名称:forcebalance,代码行数:70,代码来源:npt_tinker.py


注:本文中的forcebalance.forcefield.FF.print_map方法示例由纯净天空整理自Github/MSDocs等开源代码及文档管理平台,相关代码片段筛选自各路编程大神贡献的开源项目,源码版权归原作者所有,传播和使用请参考对应项目的License;未经允许,请勿转载。