def two_gaussian_potential_bc(vnew, f2, coords):
    """
    Apply Boundary Condition to the potential, force, and coordinates.

    Parameters:
    -----------
    vnew       : float (or array of floats)
                 Potential Energy

    f2         : float (or array of floats)
                 Force

    coords     : float
                 coordinates

    Returns:
    --------

    vnew       : float (or array of floats)
                 Adjusted potential energy from boundary condition

    F          : float (or array of floats)
                 Adjusted force from boundary condition

    coords     : float
                 adjusted coordinates from boundary condition

    bcbias     : float
                 bias applied strictly from the boundary condition
    """
    vold = vnew
    bcbias = 0
    is_periodic = False
    if (coords < -4.3193):

        vnew = 100.0 * (coords+4.0)**4.0 - 1.690133
        f2 = -100.0 * 4.0 * (coords+4.0)**3.0
        bcbias = vnew - vold

    elif (coords > 4.25882):

        vnew = 100.0 * (coords-4.0)**4.0 - 0.845067
        f2 = -100.0 * 4.0 * (coords-4.0)**3.0
        bcbias = vnew - vold

    return (vnew, f2, coords, bcbias, is_periodic)