single excitation configuration interaction ( cis),(1 ) time - dependent hartree  fock ( tdhf ) , and linear response time - dependent density functional theory ( tddft ) are widely used for ab initio calculations of electronic excited states of large molecules ( more than 50 atoms , thousands of basis functions ) because these single - reference methods are computationally efficient and straightforward to apply .
although highly correlated and/or multireference methods , such as multireference configuration interaction ( mrci),(10 ) multireference perturbation theory ( mrmp(11 ) and caspt2 ) , ( 12)and equation - of - motion coupled cluster methods ( sac - ci(13 ) and eom - cc ) , allow for more reliably accurate treatment of excited states , including those with double excitation character , these are generally too computationally demanding for large molecules .
cis / tdhf is essentially the excited - state corollary of the ground - state hartree 
fock ( hf ) method and thus similarly suffers from a lack of electron correlation . because of this
, cis / tdhf excitation energies are consistently overestimated , often by 1 ev.(8 ) the tddft method includes dynamic correlation through the exchange  correlation functional , but standard nonhybrid tddft exchange  correlation functionals generally underestimate excitation energies , particularly for rydberg and charge - transfer states.(5 ) the problem in charge - transfer excitation energies is due to the lack of the correct 1/r coulombic attraction between the separated charges of the excited electron and hole.(16 ) charge - transfer excitation energies are generally improved with hybrid functionals and also with range separated functionals that separate the exchange portion of the dft functional into long- and short - range contributions .
neither cis nor tddft ( with present - day functionals ) properly includes the effects of dispersion but promising results have been obtained with an empirical correction to standard dft functionals , and there are continued efforts to include dispersion directly in the exchange  correlation functional .
both the cis and tddft(26 ) single reference methods lack double excitations and are unable to model conical intersections or excitations in molecules that have multireference character . in spite of these limitations , the cis and tddft methods can be generally expected to reproduce trends for one - electron valence excitations , which are a majority of the transitions of photochemical interest .
tddft using hybrid density functionals , which include some percentage of hf exact exchange , has been particularly successful in modeling the optical absorption of large molecules .
furthermore , the development of new dft functionals and methods is an avid area of research , with many new functionals introduced each year .
thus it is a virtual certainty that the quality of the results available from tddft will continue to increase .
a summary of the accuracy currently available for vertical excitation energies is available in a recent study by jacquemin et al . which compares tddft results using 29 functionals for 500 molecules.(29 ) although cis and tddft are the most tractable methods for excited states of large molecules , their computational cost still prevents application to many systems of photochemical interest .
thus , there is considerable interest in extending the capabilities of cis / tddft to even larger molecules , beyond hundreds of atoms .
quantum mechanics / molecular mechanics ( qm / mm ) schemes provide a way to model the environment of a photophysically interesting molecule by treating the molecule with qm and the surrounding environment with mm force fields . however , it is difficult to know when the mm approximations break down and when a fully qm approach is necessary . with fast , large - scale cis / tddft calculations , all residues of a photoactive protein could be treated quantum mechanically to explore the origin of spectral tuning , for example .
explicit effects of solvent  chromophore interactions , including hydrogen bonding , charge transfer , and polarization , could be fully included at the ab initio level in order to model solvatochromic shifts .
one potential route to large scale cis and tddft calculations is through exploitation of the stream processing architectures(35 ) now widely available in the form of graphical processing units ( gpus ) .
the introduction of the compute unified device architecture(36 ) ( cuda ) as an extension to the c language has greatly simplified gpu programming , making it more easily accessible for scientific programming .
gpus have already been applied to achieve speed - ups of orders of magnitude in ground - state electronic structure , ab initio molecular dynamics(41 ) and empirical force field - based molecular dynamics calculations . in this article
we extend our implementation of gpu quantum chemistry in the newly developed terachem program(46 ) beyond our previous two - electron integral evaluation(47 ) and ground - state self - consistent field , geometry optimization , and dynamics calculations(41 ) to also include the calculation of excited electronic states .
we use gpus to accelerate both the matrix multiplications within the cis / tddft procedure and also the integral evaluation ( these steps comprise most of the effort in the calculation ) . the computational efficiency that arises from the use of redesigned quantum chemistry algorithms on gpu hardware to evaluate electron repulsion integrals ( eris )
allows full qm treatment of the excited states of very large systems  both large chromophores and chromophores in which the environment plays a critical role and should be treated with qm .
dancoff approximation using gpus to drastically speed up the bottleneck two - electron integral evaluation , density functional quadrature , and matrix multiplication steps .
this results in cis calculations over 400 times faster than those achieved running on a comparable cpu platform .
the linear response formalism of tdhf and tddft has been thoroughly presented in review articles .
only the equations relevant for this work are presented here , and real orbitals are assumed throughout . the tdhf / tddft working equation for determining the excitation energies  and corresponding x and y transition amplitudes iswhere for tdhf ( neglecting spin indices for simplicity):and for tddft : the two electron integrals ( eris ) are defined asand within the adiabatic approximation of density functional theory , in which the explicit time dependence of the exchange  correlation functional is neglected : the i , j and a , b indices represent occupied and virtual molecular orbitals ( mos ) , respectively , in the hf / kohn  sham ( ks ) ground - state determinant . setting the b matrix to zero within tdhf results in the cis equation , while in tddft this same neglect yields the equation known as the tamm  dancoff approximation ( tda):in part because dft virtual orbitals provide a better starting approximation to the excited state than hf orbitals , the tda generally gives results that are very close to the full linear response tddft results for nonhybrid dft functionals at equilibrium geometries .
furthermore , previous work has shown that a large contribution from the b matrix in tddft ( and to a lesser extent also in tdhf ) is often indicative of a poor description of the ground state , either due to singlet  triplet instabilities or multireference character.(52 ) casida and co - workers have examined the breakdown of tddft in modeling photochemical pathways(52 ) and have come to the conclusion that  the tda gives better results than does conventional tddft when it comes to excited - state potential energy surfaces in situations where bond breaking occurs .  thus , if there is substantial deviation between the full tddft and tda - tddft excitation energies , then the tda results will often be more accurate .
a standard iterative davidson algorithm(53 ) has been implemented to solve the cis / tda - tddft equations . as each ax matrix  vector product is formed , the required two - electron integrals are calculated over primitive basis functions within the atomic orbital ( ao ) basis directly on the gpu . within cis , the ax matrix  vector product is calculated ashere greek indices represent ao basis functions , ci is the ground - state mo coefficient of the hf / ks determinant , and t is a nonsymmetric transition density matrix .
for very small matrices , there is no time savings with gpu computation of the matrix multiplication steps in eqs 10 and 12 . for matrices of dimension less than 300
the linear algebra is performed on the gpu using calls to the nvidia cublas library , a cuda implementation of the blas library.(54 ) in quantum chemistry the ao basis functions are generally a linear combination of primitive atom - centered gaussian basis functions . for a linear combination of m primitive basis functions
 centered at a nucleus , the contracted ao basis function  isthus the two - electron integrals in the contracted ao basis that need to be evaluated for eq 11 above are given bywhere parentheses indicate integrals over contracted basis functions and square brackets indicate integrals over primitive functions .
while transfer of matrix multiplication from the cpu to the gpu provides some speedup , the gpu acceleration of the computation of the eris delivers a much more significant reduction in computer time .
details of our gpu algorithms for two - electron integrals in the j and k matrices ( coulomb and exchange operators , respectively ) have been previously published , so we only briefly highlight information relevant to our excited - state implementation , which uses these algorithms .
the gpu evaluates the j and k matrices over primitives , and these are contracted on the cpu .
initially pairs of primitive atomic orbital functions are combined using the gaussian product rule into a set of bra- and ket- pairs .
a prescreening threshold is used to remove negligible pairs , and the remaining pairs are sorted by angular momentum class and by their [ bra| or |ket ] contribution to the total [ bra|ket ] schwarz bounds , respectively.(55 ) all data needed to calculate the [ bra|ket ] integrals ( e.g. , exponents , contraction coefficients , atomic coordinates , etc . ) are then distributed among the gpus .
the coulomb j matrix and exchange k matrix are calculated separately , with both algorithms designed to minimize interthread communication by ensuring that each gpu has all necessary data for its share of integrals .
the [ bra| and |ket ] pairs are processed in order of decreasing bound , and execution is terminated once the combined [ bra|ket ] bound falls below a predetermined threshold . because the ground - state density matrix is symmetric , both the ground - state j and k matrices are also symmetric , and thus only the upper triangle of each needs to be computed .
the coulomb j matrix elements are given bywithin our j matrix algorithm , one gpu thread evaluates one primitive two - electron integral using the hermite gaussian formulation as in the mcmurchie 
davidson algorithm , which then must be contracted into the final j matrix element as given in eq 15 .
j matrix computation uses the    and    symmetry and eliminates duplicates within the bra and ket primitive hermite product lists , reducing the number of integrals calculated from o(n ) to o(n/4 ) . a different gpu subroutine ( or gpu kernel )
is called for each angular momentum class , leading to nine gpu kernel calls for all s- and p- combinations : [ ss|ss ] , [ ss|sp ] , [ ss|pp ] , [ sp|ss ] , [ sp|sp ] , [ sp|pp ] , [ pp|ss ] , [ pp|sp ] , and [ pp|pp ] .
many integrals are calculated twice because [ bra|ket ]  [ ket|bra ] symmetry is not taken into account .
this is intentional  it is often faster to carry out more ( but simpler ) computations on the gpu ( compared to an algorithm that minimizes the number of floating point operations ) in order to avoid bookkeeping overhead and/or memory access bottlenecks .
this may be viewed as a continuation of a trend that began already on cpus and has been discussed in that context previously.(58 ) the maximum density matrix element of all angular momentum components weights the ket contribution to the schwarz upper bound .
this allows the jmatrix algorithm to take complete advantage of sparsity in the density matrix , since there is a one - to - one mapping between ket pairs and density matrix elements .
also , because density matrix elements are packed together with the j matrix ket integral data , its memory access pattern is contiguous , i.e. , neighboring threads access neighboring memory addresses . in general ,
the exchange kmatrix elements are given bywithin our k matrix algorithm , one block of gpu threads evaluates one k matrix element and thereby avoids any communication with other thread blocks . because the integrals ( bra| ) and ( bra| ) are paired with different density matrix elements , the k matrix algorithm does not take into account the    and    symmetry . on the other hand , [ bra|ket ]
 [ ket|bra ] symmetry is used , leading to o(n/2 ) integrals computed to form the final k matrix .
in addition to having to compute more integrals than is required for the j matrix computation , the k matrix computation is slowed relative to j matrix computation by two additional issues .
the first is that unlike the j matrix gpu implementation , the k matrix algorithm can not map the density matrix elements onto the ket integral data , since the density index now spans both bra and ket indices .
the second issue facing k matrix computation is that because the sparsity of the density can not be included in the presorting of ket pairs , the sorted integral bounds can not be guaranteed to be strictly decreasing , and a more stringent cutoff threshold ( still based on the product of the density matrix element and the schwarz upper bound ) must be applied for k kernels , meaning that k computation does not take as much advantage of density matrix sparsity as j computation . as a result of these drawbacks , the exchange matrix takes longer to calculate than its coulomb counterpart .
based solely on the number of integrals required , the k / j timing ratio for ground - state scf calculations should be 2 . in practice , with the memory access and the thresholding issues , values of 35 are more common . in our excited - state calculations
, we use the same j and k matrix gpu algorithms , adjusted for the fact that the nonsymmetric transition density matrix t replaces the symmetric ground - state density matrix p. the portion of the f matrix from the product of t with the first integral in eq 11 is computed with the j matrix algorithm .
the portion of the f matrix from the product of t with the second integral in eq 11 is computed with the k matrix algorithm . while the j matrix remains symmetric
we must thus calculate both the upper and lower triangle contributions for the cis / tddft k matrix , resulting in two calls to the k matrix algorithm and computation of up to o(n ) integrals .
in addition to an increased number of integrals in the excited state , the k / j timing discrepancy ( comparing cis / tddft to ground - state scf calculations ) is also increased due to the sparseness of the transition density compared to the ground - state density .
evaluation of the exchange  correlation functional contribution from eq 7 needed for tddft excited states(7 ) is performed using numerical quadrature on a three - dimensional grid , which maps efficiently onto massively parallel architectures , such as the gpu .
this was recently demonstrated for ground - state dft , for both gpu and related(59 ) architectures .
the expensive steps are evaluating the electron density / gradient at the grid quadrature points to numerically determine the necessary functional derivatives and summing the values on the grid to assemble the matrix elements of eq 7 .
we use a becke - type quadrature scheme(60 ) with lebedev angular(61 ) and euler 
, we generate the second functional derivative of the exchange  correlation functional only once , saving its value at each quadrature point in memory .
then , for each davidson iteration , the appropriate integrals are evaluated , paired with the saved functional derivative values , and summed into matrix elements .
we evaluate the performance of our gpu - based cis / tddft algorithm on a variety of test systems : 6,6-bis(2-(1-triphenyl)-4-phenylquinoline ( b3ppq ) , an oligoquinoline recently synthesized and characterized by the jenekhe group for use in oled devices(63 ) and characterized theoretically by tao and tretiak;(64 ) four generations of oligothiophene dendrimers that are being studied for their interesting photophysical properties ; the entire photoactive yellow protein ( pyp)(68 ) solvated by tip3p(69 ) water molecules ; and deprotonated trans - thiophenyl - p - coumarate , an analogue of the pyp chromophore(70 ) that takes into account the covalent cysteine linkage , solvated with an increasing number of qm waters .
we use the 6 - 31 g basis set for all computations , since we do not yet have gpu integral routines implemented for d - functions .
this limits the quality of the excited - state energies , as polarization functions can give improved accuracy relative to experimental values and are often necessary for metals and hypervalent atoms , such as sulfur and phosphorus .
benchmark structures are shown in figures 1 and 2 along with the number of atoms and basis functions for a 6 - 31 g basis set . for the solvated pyp chromophore ,
only three structures are shown in figure 2 , but benchmark calculations are presented for 15 systems with increasing solvation , starting from the chromophore in vacuum and adding water molecules up to a 16  solvation shell , which corresponds to 900 water molecules .
structures , number of atoms , and basis functions ( fns ) using the 6 - 31 g basis set for four generations of oligothiophene dendrimers , s1s4 .
structures , number of atoms , and basis functions ( fns ) for the 6 - 31 g basis for benchmark systems photoactive yellow protein ( pyp ) , the solvated pyp chromophore , and oligoquinoline b3ppq . for pyp , carbon , nitrogen , oxygen , and sulfur atoms
are green , blue , red , and yellow , respectively . for the other molecules ,
atom coloration is as given in figure 1 , with additional red and blue coloration for oxygen and nitrogen atoms , respectively . for our benchmark tddft calculations ,
we use the generalized gradient approximation with becke s exchange functional(71 ) combined with the lee , yang , and parr correlation functional(72 ) ( blyp ) , as well as the hybrid b3lyp functional . during the scf procedure for the ground - state wave function , we use two different dft grids . a sparse grid of 1000 grid points / atom is used to converge the wave function until the diis error reaches a value of 0.01 , followed by a more dense grid of 3000 grid points / atom until the ground - state wave function is fully converged .
this denser grid is also used for the excited - state tddft timings reported herein , unless otherwise noted .
an integral screening threshold value of 1  10 atomic units is used by default unless otherwise noted .
within terachem , this means that coulomb integrals with products of the density element and schwarz bound below the integral screening threshold are not computed , and exchange integrals with products of the density element and schwarz bound below the threshold value times a guard factor of 0.001 are not computed .
the initial n pair quantities list is also pruned , with a default pruning value of 10 for removing pairs from integral computation .
the pair quantity pruning value is set to the smaller of 10 and 0.01  the integral screening threshold .
the timings reported herein were obtained on a desktop workstation using dual quad - core intel xeon x5570 cpus , 72 gb ram , and 8 tesla c1060 gpus .
all cpu operations are performed in full double precision arithmetic , including one - electron integral evaluation , integral postprocessing and contraction , and diagonalization of the subspace matrix of a. to minimize numerical error , integral accumulation also uses double precision .
calculations carried out on the gpu ( coulomb and exchange operator construction and dft quadrature ) use mixed precision unless otherwise noted .
the mixed precision integral evaluation is a hybrid of 32- and 64-bit arithmetic . in this case , integrals with schwarz bounds larger than 0.001 au are computed in full double precision , and all others are computed in single precision .
the potential advantages of mixed precision arithmetic in quantum chemistry have been discussed in the context of gpu architectures by several groups and stem in part from the fact that there are often fewer double precision floating point units on a gpu than single precision floating point units . to study the effects of using single precision on excited - state calculations
, we have run the same cis calculations using both single and double precision integral evaluation for many of our benchmark systems . in general
we find that mixed ( and often even single ) precision arithmetic on the gpu is more than adequate for cis / tddft . in most cases
we find that the convergence behavior is nearly identical for single and double precision until the residual vector is quite small .
figure 3 shows the typical single and double precision convergence behavior as represented by the cis residual vector norm convergence for b3ppq , the first and third generations of oligothiophene dendrimers s1 and s3 , and a snapshot of the pyp chromophore surrounded by 14 waters .
the convergence criterion of the residual norm , which is 10 au , is shown with a straight black line .
note that for the examples in figure 3 , we are not using mixed precision  all two - electron integrals on the gpu are done in single precision ( with double precision accumulation as described previously).(39 ) this is therefore an extreme example ( other calculations detailed in this paper used mixed precision where large integrals and quadrature contributions are calculated in double precision ) and serves to show that cis and tddft are generally quite robust , irrespective of the precision used in the calculation .
nevertheless , a few problematic cases have been found in which single precision integral evaluation is not adequate and where double precision is needed to achieve convergence.(75 ) during the course of hundreds of cis calculations performed on snapshots of the dynamics of the pyp chromophore solvated by various numbers of water molecules , a small number ( < 1% ) of cases yield ill - conditioned davidson convergence when single precision is used for the gpu - computed eris and quadrature contributions . for illustration , the single and double precision convergence behavior for one of these rare cases , here the pyp chromophore with 94 waters ,
is shown in figure 3 . in practice , this is not a problem since one can always switch to double precision , and this can be done automatically when convergence problems are detected .
recent work in our group(76 ) shows a speedup of 24 times for an rhf ground - state calculation in going from full double precision to mixed or single precision for our gpu eri algorithms .
plot of single and double precision ( sp and dp ) convergence behavior for the first cis/6 - 31 g excited state of five of the benchmark systems .
the convergence threshold of 10 ( norm of residual vector ) is indicated with a straight black line . in most cases ,
convergence behavior is identical for single and double precision integration until very small residual values well below the convergence threshold .
one such example is shown here for a snapshot of the pyp chromophore ( pypc ) surrounded by 94 waters .
timings and cis excitation energies ( from the ground - state s0 to the lowest singlet excited state s1 ) for some of the test systems are given in table 1 and compared to the gamess quantum chemistry package version 12 jan 2009 ( r3 ) .
the gamess timings are obtained using the same intel xeon x5570 eight - core machine as for the gpu calculations ( where gamess is running in parallel over all eight cores ) .
we compare to gamess because it is a freely available and mature quantum chemistry code and provides a reasonable benchmark of the expected speed of the algorithms on a cpu .
gamess may not represent the absolute best performance that can be achieved using the implemented algorithms on a cpu.(40 ) coordinates of all the geometries used in the tests are provided in supporting information , so the interested reader can determine timings for other codes and architectures if further comparisons are desired .
unfortunately , it is not possible to compare our own code against itself , running on the cpu or gpu , since there does not presently exist a compiler that can generate a cpu executable from cuda code .
calculations were performed on a dual intel xeon x5570 ( 8 cpu cores ) with 72 gb ram .
gpu calculations use 8 tesla c1060 gpu cards . comparing the values for the cis first excited - state energy ( e s0/s1 ) given in table 1
, we find that the numerical accuracy of the excitation energies for mixed precision gpu integral evaluation is excellent for all systems studied .
the largest discrepancy in the reported excitation energies between gamess and our gpu implementation in terachem is less than 0.00004 ev .
we also report the cis times and speedups for gamess and gpu accelerated cis in terachem ( note that the times reported refer to the entire cis calculation from the completion of the ground - state scf to the end of program execution ) .
since cis is necessarily preceded by a ground - state scf calculation , we also report the scf speedups to give a complete picture .
we leave out the absolute scf times , since the efficiency of the gpu - based scf algorithm has been discussed for other test molecules previously .
we find a large increase in performance is obtained using the gpu for both ground- and excited - state methods .
the speedups increase as system size increases , with scf speedups outperforming cis speedups . for the largest system compared with gamess , which is the 29 atom chromophore of pyp surrounded by 487 qm water molecules ,
some possible reasons for the differing speedups in ground- and excited - state calculations are discussed below . in the supporting information ,
we also include a table giving the absolute terachem scf and cis times for four of the test systems , along with the corresponding scf and cis energies , for both mixed and double precision computation and for three different integral screening threshold values .
while the timings increase considerably in switching from mixed precision to double precision and in tightening the integral screening thresholds , the cis excitation energies remain nearly identical , suggesting that the cis algorithm is quite robust with respect to thresholding .
the dominant computational parts in building the cis / tddft ax vector can be divided into coulomb j matrix , exchange k matrix , and dft contributions .
figure 4 plots the cpu + gpu time consumed by each of these three contributions ( both cpu and gpu times are included here , although the cpu time is a very small fraction of the total ) , in which j and k timings are taken from an average of the 10 initial guess ax builds for a cis calculation , and the dft timings are from an average of the initial guess ax builds for a td - blyp calculation . the initial guess
transition densities are very sparse , and thus this test highlights the differing efficiency of screening and thresholding in these three contributions .
the j timings for cis and blyp are similar , and only those for cis are reported .
power law fits are shown as solid lines and demonstrate near - linear scaling behavior of all three contributions to the ax build .
the j matrix and dft quadrature steps are closest to linear scaling , with observed scaling of n for both contributions , where n is the number of basis functions . the k matrix contribution scales as n because it is least able to exploit the sparsity of the transition density matrix .
these empirical scaling data demonstrate that with proper sorting and integral screening , the ax build in cis and tddft scales much better than quadratic , with no loss of accuracy in excitation energies .
contributions to the time for building an initial ax vector in cis and td - blyp .
ten initial x vectors are created based on the mo energy gap , and the timing reported is the average time for building ax for those 10 vectors . the timings are obtained on a dual intel xeon x5570 platform with 72 gb ram using 8 tesla c1060 gpus . data ( symbols ) are fit to power law ( solid line , fitting parameters in inset ) .
fewer points are included for the td - blyp timings because the scf procedure does not converge for the solvated pyp chromophore with a large number of waters or for the full pyp protein .
of the three integral contributions ( j , k , and dft quadrature ) , the computation of the k matrix is clearly the bottleneck .
this is due to the three issues with exchange computation previously discussed : ( 1 ) the j matrix takes full advantage of density sparsity because of efficient density screening that is not possible for our k matrix implementation , ( 2 ) exchange kernels access the density in memory noncontiguously , and ( 3 ) exchange requires the evaluation of 4 times more integrals than j both because it lacks the    and    symmetry and because it needs to be called twice to account for the nonsymmetric excited - state transition density matrix .
it is useful to compare the time required to calculate the k matrix contribution to the first ground - state scf iteration ( which is the most expensive iteration due to the use of fock matrix updating(77 ) ) and to the ax vector build for cis ( or td - b3lyp ) .
we find that for the systems studied herein the k matrix contribution is on average almost 2 times faster in cis compared to the first iteration of the ground - state scf .
one might have expected the excited - state computation to be 2 times slower because of the two k matrix calls , but the algorithm efficiently exploits the greater sparsity of the transition density matrix ( compared to the ground - state density matrix ) . due to efficient prescreening of the density and integral contributions to the schwarz bound before the gpu coulomb kernels are launched , the j matrix computation also exploits the greater sparseness of the transition density and therefore is 3.5 times faster than the ground - state first iteration j matrix computation . since j matrix computation profits more from transition density sparsity than k matrix computation , the current implementation of the j matrix computation scales better with system size than the implementation of the k matrix computation ( n vs n for the excited - state benchmarks presented here ) . as can be seen in figure 4,(78 ) the dft integration usually takes more time than the j matrix contribution .
this is because of the larger prefactor for dft integration , which is related to the density of the quadrature grids used .
it has previously been noted(79 ) that very sparse grids can be more than adequate for tddft .
we further support this claim with the data presented in table 2 , where we compare the lowest excitation energies and the average td - blyp integration times for the initial guess vectors for six different grids on two of the test systems . for both molecules ,
the excitation energies from the sparsest grid agree well with those of the more dense grids but with a substantial reduction in integration time , suggesting that a change to an ultra sparse grid for the tddft portion of the calculation could result in considerable time savings with little to no loss of accuracy .
the small ( < 0.0002 ev ) differences in excitation energies between our gpu - based td - blyp and the cpu - based nwchem are likely due to slightly differing ground - state densities , which differ in energy by 7 microhartrees for the chromophore and 1.9 millihartrees for the s2 dendrimer .
td - blyp timings ( average time for the dft quadrature in one ax build for the initial 10 ax vectors ) . for comparison ,
number of points / atom refers to the pruned grid for terachem and the unpruned grid for nwchem .
nwchem was run on a different architecture , so timings are not comparable . while successive ground - state scf iterations take less computation time than the first ( because of the use of fock matrix updating ) , all iterations in the excited - state calculations take roughly the same amount of time .
this is the dominant reason for the discrepancy in the speedups for ground - state scf and excited - state cis shown in table 1 .
an additional reason that the scf speedup is greater than the cis speedup is decreased parallel efficiency because the ground - state density is less sparse than the transition density ( all of the reported calculations are running on eight gpu cards in parallel ) .
gpu - accelerated cis and tddft computation provides the excited states of much larger compounds than can be currently studied with ab initio methods . for the well - behaved valence transitions in the pyp systems
, cis convergence requires very few davidson iterations . the total wall time ( scf + cis ) required to calculate
the first cis/6 - 31 g excited state of the entire pyp protein ( 10869 basis functions ) is less than 6 h , with 4.7 h devoted to the scf procedure and 1.2 h to the cis procedure .
we can thus treat the protein with full qm and study how mutation within pyp will affect the absorbance . for any meaningful comparison with the experimental absorption energy of pyp at 2.78 ev,(70 ) many geometrical configurations need to be taken into account . for this single configuration ,
the cis excitation energy of 3.69 ev is much higher than the experimental value , as expected with cis .
the td - b3lyp bright state ( s5 ) is closer to the experimental value but still too high at 3.33 ev .
solvatochromic studies in explicit water are problematic for standard dft methods , including hybrid functionals , due to the well - known difficulty in treating charge - transfer excitations . in calculating the timings for the first excited state of the pyp chromophore with increasing numbers of waters
, we found that the energy of the cis first excited state quickly leveled off and stabilized , while that for td - blyp and td - b3lyp generally decreased to unphysical values , at which point the ground - state scf convergence was also problematic .
this behavior of the first excitation energies for the pyp chromophore with increasing numbers of waters is shown in figure 5 for cis , td - blyp , and td - b3lyp .
while the 20% hf exchange in the hybrid td - b3lyp method does improve the excitation energies over td - blyp , the energies are clearly incorrect for both methods , and a higher level of theory or a range - separated hybrid functional is certainly necessary for studying excitations involving explicit qm waters . the first excitation energy ( ev ) of the pyp chromophore with increasing numbers of surrounding water molecules .
both td - blyp and td - b3lyp exhibit spurious low - lying charge - transfer states .
the recent theoretical work by badaeva et al . examining the one and two photon absorbance of oligothiophene dendrimers was limited to results for the first three generations s1s3 , even though experimental results were available for s4 . in table 3
, we compare our gpu accelerated results on the first bright excited state ( oscillator strength > 1.0 ) using td - b3lyp within the tda to the full td - b3lyp and experimental results .
results within the tda are comparable to those from full td - b3lyp , for both energies and transition dipole moments . our results for s4 show the continuing trend of decreasing excitation energy and increasing transition dipole moment with increasing dendrimer generation .
we have implemented ab initio cis and tddft calculations within the terachem software package , designed from inception for execution on gpus .
the numerical accuracy of the excitation energies is shown to be excellent using mixed precision integral evaluation .
the ability to use lower precision in much of the cis and tddft calculation is reminiscent of the ability to use coarse grids when calculating correlation energies , as shown previously for pseudospectral methods .
recently , it has also been shown(86 ) that single precision can be adequate for computing correlation energies with cholesky decomposition methods which are closely related to pseudospectral methods.(87 ) both quadrature and precision errors generally behave as relative errors , while chemical accuracy is an absolute standard ( often taken to be 1 kcal / mol ) . thus , coarser grids and/or lower precision can be safely used when the quantity being evaluated is itself small ( and therefore less relative accuracy is required ) , as is the case for correlation and/or excitation energies . for some of the smaller benchmark systems , we present speedups as compared to the gamess quantum chemistry package running over eight processor cores .
the speedups obtained for cis calculations range from 9 to 461 times , with increasing speedups with increasing system size .
these speedup figures are not necessarily normative ( other quantum chemistry packages might be more efficient ) , but we feel they give a good sense of the degree to which redesign of quantum algorithms for gpus may be useful . the increased size of the molecules that can be treated using our gpu - based algorithms exposes some failings of dft and tddft .
specifically , the charge - transfer problem(16 ) of tddft and the delocalization problem(88 ) of dft both seem to become more severe as the molecules become larger , especially for the case of hydrated chromophores with large numbers of surrounding quantum mechanical water molecules .
it remains to be seen whether range - separated hybrid functionals can solve these problems for large molecules , and we are currently working to implement these .
 excited - state calculations are implemented in a development version of the gpu - based terachem software package using the configuration interaction singles ( cis ) and adiabatic linear response tamm  
 dancoff time - dependent density functional theory ( tda - tddft ) methods . 
 the speedup of the cis and tddft methods using gpu - based electron repulsion integrals and density functional quadrature integration allows full ab initio excited - state calculations on molecules of unprecedented size . 
 cis/6 - 31 g and td - blyp/6 - 31 g benchmark timings are presented for a range of systems , including four generations of oligothiophene dendrimers , photoactive yellow protein ( pyp ) , and the pyp chromophore solvated with 900 quantum mechanical water molecules . 
 the effects of double and single precision integration are discussed , and mixed precision gpu integration is shown to give extremely good numerical accuracy for both cis and tddft excitation energies ( excitation energies within 0.0005 ev of extended double precision cpu results ) . 
Introduction
CIS/TDDFT Implementation using GPUs
Results and Discussion
Conclusions