A new robust method for variational determination of atomic zero-flux surfaces is presented. The zero-flux surface sheets are expressed in terms of variational trial. functions in prolate spheroidal coordinates. The trial functions are optimized with a Newton procedure to satisfy the zero-flux condition on a grid. The data required for radial integrations are generated by an adaptive quadrature procedure that employs model electron densities and utilizes an original third-order algorithm for Linear search. Results of test calculations involving variational determination of atomic surfaces are presented for a representative set of 20 molecules. The new approach is both less time consuming and substantially more accurate than the previously published algorithms. (C) 1995 by John Wiley & Sons, Inc.