We present a novel parallel gradient optimization algorithm designed for the optimization of molecular geometry - the parallel preconditioned LBFGS (PP-LBFGS) method. In each step, several additional gradient calculations (performed in parallel with the calculation of the potential) are used to improve the most important elements of the Hessian. The sparsity of the connectivity matrix and the graph theory are used to estimate multiple Hessian elements from each additional gradient calculation. The simplest variant of the algorithm, which requires 4 gradient evaluations per cycle, converges 2x-4x faster than the LBFGS algorithm, depending on the size of the system. (C) 2013 Elsevier B. V. All rights reserved.