import math


def torsion_energy(phi, fc, multiplicity, npaths, phi0):
    """
    Calculate the torsion energy using the periodic potential.

    One torsion term at a time, using this function:
    E = fc/npaths * [1 + cos(multiplicity*phi - phi0)]

    Args:
      phi (float):  dihedral angle [degrees]
      fc (float):  force constant [kcal/mol]
      multiplicity (float):  multiplicity [/]
      npaths (float):  division factor [/]
      phi0 (float):  phase shift [degrees]

    Returns:
      e_torsion (float):  energy of torsion [kcal/mol]
    """
    phi = math.pi/180.0 * phi   # degrees to radians
    phi0 = math.pi/180.0 * phi0   # degrees to radians
    return 1.0 * fc/npaths * (1 + math.cos(multiplicity*phi - phi0))