The COLUMBUS program system is a collection of Fortran programs for performing general multireference single‐ and double‐excitation configuration interaction (MRSDCI) wave function optimization based on the graphical unitary group approach. The program system also includes integral generation, SCF and MCSCF orbital optimization, integral transformation, and wave function analysis programs. The original program system was written in 1980 to 1981. Since that time, it has evolved into one of the most popular MRSDCI program systems used in the computational chemistry community. The discussion of this evolution will include the exploitation of efficient matrix‐matrix and matrix‐vector product computational kernels, the use of generally contracted symmetry‐adapted orbital basis sets, general Hamiltonian diagonalization procedures, energy‐based internal walk selection, flexible DRT specification, improved coupling‐coefficient evaluation methods, coupled‐pair functional and multireference CPF capabilities, and density matrix construction. The numerous versions of the program system that are maintained at different sites and on different computers are now in the process of being merged. The source code for this combined version will be made available to the computational chemistry community. The source code for a specific computer may be generated from the source code for another computer by a single pass through a simple filter utility that is included with the program system. The directly supported computers will initially include various models of VAX, Cray, FPS, IBM, CDC, and ETA machines with the addition of other machines shortly thereafter. The ongoing developments of the COLUMBUS system that are discussed include a new method for computing analytic energy gradients for MRSDCI wave functions. This effective‐density‐matrix based method avoids the “coupled perturbed MCSCF” solutions for each coordinate direction, avoids the transformation of any derivative‐integral quantities from the AO to the MO basis, avoids the transformation of the coupling coefficients from the MO to the AO basis, allows a subset of the MCSCF doubly occupied orbitals to be frozen in the CI wave function, and allows the MRSDCI wave function to be generated from general reference CSFs that are not necessarily related to the MCSCF expansion CSFs.