A Parallel Algorithm for Solving the 3d Schrödinger Equation
Abstract
We describe a parallel algorithm for solving the timeindependent 3d Schrödinger equation using the finite difference time domain (FDTD) method. We introduce an optimized parallelization scheme that reduces communication overhead between computational nodes. We demonstrate that the compute time, , scales inversely with the number of computational nodes as . This makes it possible to solve the 3d Schrödinger equation on extremely large spatial lattices using a small computing cluster. In addition, we present a new method for precisely determining the energy eigenvalues and wavefunctions of quantum states based on a symmetry constraint on the FDTD initial condition. Finally, we discuss the usage of multiresolution techniques in order to speed up convergence on extremely large lattices.
pacs:
03.65.Ge, 02.70.c, 02.30.Jr, 02.70.BfI Introduction
Solving the 3d Schrödinger equation given an arbitrary potential is of great practical use in modern quantum physics; however, there are only a handful of potentials for which analytic solution is possible. In addition, any potential that does not have a high degree of symmetry, e.g. radial symmetry, requires solution in full 3d, making standard “pointandshoot” methods numrepc () for solving onedimensional partial differential equations of little use. In this paper we discuss a parallel algorithm for solving the 3d Schrödinger equation given an arbitrary potential using the finite difference time domain (FDTD) method.
The FDTD method has a long history of application to computational electromagnetics Yee:1966 (); Taflove:1980 (); Taflove:1990 (); Taflove:2005 (). In the area of computational electromagnetics parallel versions of the algorithms have been developed and tested Varadarajan:1994 (); Liu:1994 (); Tinniswodd:1996 (); Sypniwski:2000 (); Palaniappan:2000 (); Yu:2005 (); Yu:2007 (). In this paper, we discuss the application of parallelized FDTD to the 3d Schrödinger equation. The standard FDTD method has been applied to the 3d Schrödinger equation by several authors in the past feagin (); roy1 (); wadehra (); sullivan1 (); sullivan2 (); roy2 (); Sudiarta:2007 (). Here we show how to efficiently parallelize the algorithm. We describe our parallel algorithm for finding ground and excited state wavefunctions and observables such as energy eigenvalues, and rootmeansquared radii. Additionally, we introduce a way to use symmetry constraints for determining excited state wavefunctions/energies and introduce a multiresolution technique that dramatically decreases compute time on large lattices. This paper is accompanied by an opensource release of a code that implements the algorithm detailed in this paper. The code uses the Message Passing Interface (MPI) protocol for message passing between computational nodes.
We note that another popular method for numerical solution of the 3d Schrödinger equation is the Diffusion Monte Carlo (DMC) technique, see Grimm:1971 (); Anderson:1975 (); Lester:1990 (); Mitas:1996 (); Foulkes:2001 () and references therein. The starting point for this method is the same as the FDTD method applied here, namely transformation of the Schrödinger equation to imaginary time. However, in the DMC algorithm the resulting “dynamical” equations are transformed into an integral Green’s function form and then the resulting integral equation is computed using stochastic sampling. The method is highly inefficient unless importance sampling Metropolis:1953 (); Umrigar:1993 () is used. DMC is efficiently parallelized and there are several codes which implement parallelized DMC QMCBeaver:2004 (); QUMAX:2004 (); AspuruGuzik:2005 (). The method is similar in many ways to the one presented herein; however, the method we use does not suffer from the fermion sign problem which forces DMC to use the socalled “fixednode approximation” Anderson:1975 (). In addition, although the DMC algorithm can, in principle, be applied to extract properties of the excited states of the system most applications to date only calculate the ground state wavefunction and its associated expectation values. The FDTD method described herein can extract both ground and excited state wavefunctions.
The organization of the paper is as follows. In Secs. II and III we briefly review the basics of the FDTD method applied to the 3d Schrödinger equation and derive the equations necessary to evolve the quantummechanical wavefunction. In Sec. IV we discuss the possibility of imposing a symmetry constraint on the FDTD initial condition in order to pick out different quantummechanical states. In Sec. V we describe our strategy for parallelizing the FDTD evolution equations and the measurement of observables. In Sec. VI we introduce an efficient method of using lowerresolution FDTD wavefunctions as initial conditions for higherresolution FDTD runs that greatly speeds up determination of highaccuracy wavefunctions and their associated observables. In Sec. VII we give results for a few potentials including benchmarks showing how the code scales as the number of computational nodes is increased. Finally, in Sec. VIII we conclude and give an outlook for future work.
Ii Setup and Theory
In this section we introduce the theory necessary to understand the FDTD approach for solving the timeindependent Schrödinger equation. Here we will briefly review the basic idea of the FDTD method and in the next section we will describe how to obtain the discretized “equations of motion”.
We are interested in solving the timeindependent Schrödinger equation with a static potential and a particle of mass
(1) 
where is a quantummechanical wavefunction that solves this equation, is the energy eigenvalue corresponding to , and is the Hamiltonian operator. In order to solve this timeindependent (static) problem it is efficacious to consider the timedependent Schrödinger equation
(2) 
A solution to (2) can be expanded in terms of the basis functions of the timeindependent problem, i.e.
(3) 
where are expansion coefficients which are fixed by initial conditions ( represents the ground state, the first excited state, etc.) and is the energy associated with each state.^{1}^{1}1The index is understood to represent the full set of quantum numbers of a given state of energy . In the degenerate case is an admixture of the different degenerate states.
By performing a Wick rotation to imaginary time, , and setting and in order to simplify the notation, we can rewrite Eq. (2) as
(4) 
which has a general solution of the form
(5) 
Since , for large imaginary time the wavefunction will be dominated by the ground state wavefunction . In the limit goes to infinity we have
(6) 
Therefore, if one evolves Eq. (4) to large imaginary times one will obtain a good approximation to the ground state wavefunction.^{2}^{2}2In this context a large imaginary time is defined relative to the energy splitting between the ground state and the first excited state, e.g. ; therefore, one must evolve to imaginary times much larger than .
This allows one to determine the ground state energy by numerically solving equation (4) for large imaginary time, and then use this wavefunction to find the energy expectation value :
(7) 
However, the method is not limited to extraction of only the ground state wavefunction and expectation values. In the next sections we will describe two different methods that can be used to extract, in addition, excited state wavefunctions.
Iii The Finite Difference Time Domain Method
To numerically solve the Wickrotated Schrödinger equation (4) one can approximate the derivatives by using discrete finite differences. For the application at hand we can, without loss of generality, assume that the wavefunction is realvalued as long as the potential is realvalued. The imaginary time derivative becomes
(8) 
where is some finite change in imaginary time.
Similarly, the right hand side of equation (4) becomes
(9)  
where, in the last term, we have averaged the wavefunction in imaginary time in order to improve the stability of the algorithm following Taflove Taflove:1990 () and Sudiarta and Geldart Sudiarta:2008 (). Note that if the potential has singular points these have to be regulated in some way, e.g. by ensuring that none of the lattice points coincides with a singular point. Assuming, for simplicity, that the lattice spacing in each direction is the same so that this equation can be rewritten more compactly by defining a difference vector
(10) 
together with a matrixvalued field
(11) 
giving
(12)  
Rewriting equation (4) with equations (8) and (12) gives the following update equation for in imaginary time:
(13) 
where and are
(14) 
and we have reintroduced the mass, , for generality. Evolution begins by choosing a random 3d wavefunction as the initial condition. In practice, we use Gaussian distributed random numbers with an amplitude of one. The boundary values of the wavefunction are set to zero; however, other boundary conditions are easily implemented. ^{3}^{3}3We note that parallel implementations of absorbing boundary conditions may present a bottleneck for the parallel calculation. Many implementations of perfectly matched layers exist Berenger:1994 (); Chew:1994 (); Gedney:1996 (); however, only a few efficient parallel implementations exist with the maximum efficiency of achieved using the WEPML scheme Zhou:2001 (); Rickhard:2002 (); Ramadan:2003 (); Ramadan:2007 (); Ramadan:2008 (). To the best of our knowledge the PML method has only been applied to unparallelized solution of the Schrödinger equation Fevens:1999 (); AlonsoMallo:2003 (); Han:2004 (). Note that during the imaginary time evolution the norm of the wavefunction decreases (see Eq. 5), so we additionally renormalize the wavefunction during the evolution in order to avoid numerical underflow. This does not affect physical observables.
We solve Eq. (13) on a threedimensional lattice with lattice spacing and lattice sites in each direction. Note that the lattice spacing and size should be chosen so that the states one is trying to determine (i) fit inside of the lattice volume, i.e. , and (ii) are described with a sufficiently fine resolution, i.e. . Also note that since we use an explicit method for solving the resulting partial differential equation for the wavefunction, the numerical evolution in imaginary time is subject to numerical instability if the time step is taken too large. Performing the standard von Neumann stability analysis Mitchell:1980 () one finds that in order achieve stability. For a fixed lattice volume , therefore, when keeping the lattice volume fixed. The total compute time scales as and assuming , we find that the total compute time scales as .
At any imaginary time the energy of the state, , can be computed via a discretized form of equation (7)
(15) 
Excited states are extracted by saving the full 3d wavefunction to local memory periodically, which we will call taking a “snapshot” of the wavefunction. After convergence of the ground state wavefunction these snapshots can be used, one by one, to extract states with higherenergy eigenvalues by projecting out the ground state wavefunction, then the first excited state wavefunction, and so on Sudiarta:2007 (). In principle, one can extract as many states as the number of snapshots of the wavefunction saved during the evolution. For example, assume that we have converged to the ground state and that we also have a snapshot version of the wavefunction taken during the evolution. To extract the first excited state we can project out the ground state using
(16) 
For this operation to give a reliable approximation to the snapshot time should obey . One can use another snapshot wavefunction that was saved and obtain the second excited state by projecting out both the ground state and the first excited state.
Finally we mention that one can extract the binding energy of a state by computing its energy and subtracting the value of the potential at infinity
(17) 
where with as usual. Note that if is a constant, then Eq. (17) simplifies to .
Iv Imposing symmetry conditions on the initial wavefunction
Another way to calculate the energies of the excited states is to impose a symmetry constraint on the initial conditions used for the FDTD evolution. The standard evolution calls for a random initial wavefunction; however, if we are solving a problem that has a potential with sufficient symmetry we can impose a symmetry condition on the wavefunction in order to pick out the different states required. For example, if we were considering a spherically symmetric Coulomb potential then we could select only the , , , etc. states by requiring the initial condition to be reflection symmetric about the , , and axes.^{1}^{1}1Of course, for a spherically symmetric potential a fully 3d method for solving the Schrödinger equation is unecessary since one can reduce the problem to solving a 1d partial differential equation. We only use this example because of its familiarity. This would preclude the algorithm finding any antisymmetric states such as the state since evolution under the Hamiltonian operator cannot break the symmetry of the wavefunction. Likewise to directly determine the excited state one can start by making the FDTD initial state wavefunction antisymmetric about one of the axes, e.g. the axis. As we will show below this provides for a fast and accurate method for determining the lowlying excited states.
Notationally, we will introduce two symbols, the symmetrization operator and the antisymmetrization operator . Here labels the spatial direction about which we are (anti)symmetrizing, i.e. . Although not required, it is implicit that we perform the symmetrization about a plane with , , or , respectively. In practice these are implemented by initializing the lattice and then simply copying, or copying plus flipping the sign, elements from one half of the lattice to the other. In practice, we find that due to roundoff error one should reimpose the symmetry condition periodically in order to guarantee that lowerenergy eigenstates do not reappear during the evolution.
V FDTD Parallelization Strategy
Parallelizing the FDTD algorithm described above is relatively straightforward. Ideally, one would segment the volume into equal subvolumes and distribute them equally across all computational nodes; however, in this paper we will assume a somewhat simpler possibility of dividing the lattice into “slices”. Our method here will be to start with a lattice and slice it along one direction in space, e.g. the direction, into pieces where is divisible by . We then send each slice of lattice to a separate computational node and have each computational node communicate boundary information between nodes which are evolving the sublattices to its right and/or left. The partitioning of the lattice is indicated via a 2d sketch in Fig. 1. In practice, in order to implement boundary conditions and synchronization of boundaries between computation nodes compactly in the code, we add “padding elements” to the overall lattice so that the actual lattice size is . The outside elements of the physical lattice hold the boundary value for the wavefunction. In all examples below the boundary value of the wavefunction will be assumed to be zero; however, different types of boundary conditions are easily accomodated. When slicing the lattice in order to distribute the job to multiple computational nodes we keep padding elements on each slice so that the actual size of the slices is . Padding elements on nodes that share a boundary are used to keep them synchronized, while padding elements on nodes that are at the edges of the lattice hold the wavefunction boundary condition.
In Fig. 2 we show a flow chart that outlines the basic method we use to evolve each node’s sublattice in imaginary time. In the figure each column corresponds to a separate computational node. Solid lines indicate the process flow between tasks and dashed lines indicate data flow between computational nodes. Shaded boxes indicate nonblocking communications calls that allow the process flow to continue while communications take place. As can be seen from Fig. 2 we have optimized each lattice update by making the first step in each update iteration a nonblocking send/receive between nodes. While this send/receive is happening each node can then update the interior of its sublattice. For example, in the two node case show in Fig. 1 this means that node 1 would update all sites with an index between 1 and 3 while node 2 would update sites with index between 6 and 8. Once these interior updates are complete each node then waits for the boundary communication initiated previously to complete, if it has not already done so. Once the boundaries have been synchronized, the boundary elements themselves can be updated. Going back to our example shown in Fig. 1 this would mean that node 1 would update all sites with index of 4 and node 2 would update all sites with an index of 5.
Convergence is determined by checking the ground state binding energy periodically, e.g. every one hundred time steps, to see if it has changed by more than a given tolerance. In the code, the frequency of this check is an adjustable parameter and should be tuned based on the expected energy of the state, e.g. if the energy is very close to zero then convergence can proceed very slowly and the check frequency should be correspondingly larger. Parametrically the check frequency should scale as .
For computation of observables each computational node computes its contribution to the observable. Then a parallel call is placed that collects the local values computed into a central value stored in computational node 1. Then node 1 broadcasts the value to the other nodes so that all nodes are then aware of the value of the particular observable. For example, to compute the energy of the state as indicated in Eq. (15) each computational node computes the portion of the sum corresponding to its sublattice and then these values are collected via a parallel sum operation to node 1 and then broadcast out to each node. Each node can then use this information to determine if the wavefunction evolution is complete. We note that the normalization of the wavefunction is done in a similar way with each node computing its piece of the norm, collecting the total norm to node 1, broadcasting the total norm to all nodes, and then each node normalizes the values contained on its sublattice. In this way computation of observables and wavefunction normalization is also parallelized in our approach.
v.1 Scaling of our 1d partitioning
In most settings computational clusters are limited by their communication speed rather than by CPU speed. In order to understand how things scale we introduce two time scales: which is the amount of time needed to update one lattice site and which is the amount of time needed to communicate (send and receive) the information contained on one lattice site. Typically unless the cluster being used is on an extremely fast network. Therefore, the algorithm should be optimized to reduce the amount of communications required.
For the onedimensional partitions employed here
(18) 
where is the number of 1d slices distibuted across the cluster and the factor of 2 comes from the 2 surfaces which must be communicated by the internal partitions. For the calculation not to have a communications bottleneck we should have . Using (18) we find that this constraint requires
(19) 
In the benchmarks section below we will present measurements of and using our test cluster. We find that . Using this, and assuming, as a concrete example, a lattice size of we find . For clusters with more than 102 nodes it would be more efficient to perform a fully 3d partitioning. In the case of a fully 3d partitioning one finds that the limit due to communications overhead is .
Vi The MultiResolution Technique
If one is interested in highprecision wavefunctions for lowlying excited states, an efficient way to do this is to use a multiresolution technique. This simply means that we start with a random wavefunction on small lattice, e.g. , and use the FDTD technique to determine the ground state and first few excited states and save the wavefunctions, either in local memory or disk. We can then use a linear combination of the coarse versions of each state as the initial condition on a larger lattice, e.g. , while keeping the lattice volume fixed. We can then “bootstrap” our way up to extremely large lattices, e.g. on the order of , by proceeding from low resolution to high resolution. In the results section we will present quantitative measurements of the speed improvement that is realized using this technique.
Vii Results
In this section we present results obtained for various 3d potentials and benchmarks that show how the code scales with the number of computational nodes. Our benchmarks were performed on a small cluster of 4 servers, each with two quadcore 2 GHz AMD Opteron processors. Each server can therefore efficiently run eight computational processes simultaneously, allowing a maximum of 32 computational nodes. ^{4}^{4}4Due to a server upgrade during publishing we were able to extend to 64 computational nodes in the general benchmark section. The servers were networked with commercial 1 Gbit/s TCP/IP networking. For the operating system we used 64 bit Ubuntu Server Edition 8.10 Linux.
vii.1 Implementation
In order to implement the parallel algorithm we use a mixture of C/C++ and the Message Passing Interface (MPI) library for message passing between computational nodes mpi (). The servers used the OpenMPI implementation of the MPI API.^{5}^{5}5The code was also tested against the MPICH implementation of the MPI API with similar results. The code itself is opensourced under the Gnu General Public License (GPL) and is available for internet download via the URL in Ref. code ().
vii.2 General Benchmarks
In this section we present data for the scaling of the time of one iteration and the time for communication on a lattice. As discussed in Sec. V.1 we expect to see ideal scaling of the code as long as communication time is shorter than the update time, i.e. . In Fig. 3a we show the time to complete one iteration as a function of the number of computational nodes on a loglog axis along with a linear fit. The linear fit obtained gives . In addition, in Fig. 3b we show a comparison of the full time for each iteration with the amount of time needed to communicate a lattice site’s information (in this case the local value of the wavefunction). In both Fig. 3a and 3b the error bars are the standard error determined by averaging over 10 runs, , where is the standard deviation across the sampled set of runs and is the number of runs.
As can be seen from Fig. 3b using a lattice the algorithm performs well up to at which point the communication time becomes equal to the iteration time. For we would see a violation of the scaling above due to communication overhead. Note that this is rough agreement with our estimate from Sec. V.1 which, for lattice predicts the point where communications and update times to be equal to be . Note that in Fig. 3b the increase in communication times as increases is due to the architecture of the cluster used for the benchmarks which has eight cores per server. If then all jobs run on one server, thereby decreasing the communications overhead. In the next section, we will present benchmarks for different potentials in order to (a) confirm the scaling obtained above in specific cases and (b) to verify that the code converges to the physically expected values for cases which are analytically solvable.
vii.3 Coulomb Potential Benchmarks
We use the following potential for finding the Coulomb wavefunctions
(20) 
where is the lattice spacing in units of the Bohr radius and is the distance from the center of the 3d lattice. The constant of is added for in order to ensure that the potential is continuous at . This is equivalent to making the potential constant for and shifting the entire potential by a constant which does not affect the binding energy. Analytically, in our natural units the binding energy of the th state is where is the principal quantum number labeling each state. The ground state therefore has a binding energy of and the first excited state has , etc. Note that to convert these to electron volts you should multiply by 27.2 eV.
In Fig. 4 we show the amount of time needed in seconds to achieve convergence of the ground state binding energy to a part in as a function of the number of computational nodes for on a loglog plot. For this benchmark we used a lattice with , a constant lattice spacing of , a constant imaginary time step of , and the particle mass was also set to . In order to remove runbyrun fluctuations due to the random initial conditions we used the same initial condition in all cases. In Fig. 4 the error bars are the standard error determined by averaging over 10 runs, , where is the standard deviation across the sampled set of runs and is the number of runs. In all cases shown the first two energy levels obtained were and . This corresponds to an accuracy of 0.2% and 2.4%, respectively. In Fig. 4 the extracted scaling slope is close to 1 indicating that the compute time in this case scales almost ideally, i.e. inversely proportional to the number of computing nodes. Note that the fit obtained in Fig. 4 has a slope with magnitude greater than 1 indicating scaling which is better than ideal; however, as one can see from the figure there is some uncertainty associated with this fit.
vii.4 3d Harmonic Oscillator Benchmarks
We use the following potential for finding the 3d harmonic oscillator wavefunctions
(21) 
where is the distance from the center of the 3d lattice.
In Fig. 5 we show the amount of time needed in seconds to achieve convergence of the ground state binding energy to a part in as a function of the number of computational nodes for . For this benchmark we used a constant lattice spacing of , a constant imaginary time step of , and a dimension lattice so that the box dimension was . In Fig. 5 the error bars are the standard error determined by averaging over 10 runs, , where is the standard deviation across the sampled set of runs and is the number of runs. The particle mass was also set to . In order to remove runbyrun fluctuations due to the random initial conditions we used the same initial condition in all cases. In all cases the ground state energy obtained was corresponding to an accuracy of 0.0026%. In Fig. 5 the extracted scaling slope is 0.91 meaning that the compute time scales as in this case. This is a slightly different slope than in the Coulomb potential case. This is due to fluctuations in compute time due to server load and sporadic network delays. The scaling coefficient reported in the conclusions will be the average of all scaling coefficients extracted from the different potentials detailed in this paper.
vii.5 Dodecahedron Potential
The previous two examples have spherical symmetry and hence it is not necessary to apply a fully 3d Schrödinger equation solver to them. We do so only in order to show scaling with computational nodes and percent error compared to analytically available solutions. As a nontrivial example of the broad applicability of the FDTD technique we apply it to a potential that is a constant negative value of inside a surface defined by a regular dodecahedron with the following 20 vertices
(22) 
where is the golden ratio. The value 1 is mapped to the point and the value 1 is mapped to the point in all three dimensions. As a result, the containing sphere has a radius of .
In Fig. 6 we show the ground and first excited states extracted from a run on a lattice with a lattice spacing of , an imaginary time step of and particle mass of . On the left we show the ground state and on the right the first excited state. We find that the energies of these two levels are and . Note that for the first excited state the position of the node surface can change during each run due to the random initial conditions used. In practice, the node surface seems to align along one randomly chosen edge of one of the pentagons that make up the surface of the dodecahedron.
In Fig. 7 we show the amount of time needed in seconds to achieve convergence of the dodecahedron ground state binding energy to a part in as a function of the number of computational nodes for . For this benchmark we used a lattice with a constant lattice spacing of , an imaginary time step of and particle mass of . In Fig. 7 the error bars are the standard error determined by averaging over 10 runs, , where is the standard deviation across the sampled set of runs and is the number of runs. In all cases the ground state energy obtained was . In Fig. 7 the extracted scaling slope is 0.91 meaning that the compute time scales as in this case.
vii.6 Applying Symmetry Constraints to the FDTD initial wavefunction
One of the fundamental problems associated with using a single FDTD run to determine both the ground state and excited states is that typically the excited states are much more extended in space than the ground state, particularly for potentials with a “long range tail” like the Coulomb potential. For this reason it is usually difficult to obtain accurate energy eigenvalues for both ground and excited states unless the lattice has an extremely fine lattice spacing and a large number of points in each direction so that the dimension of the box is also large. In Sec. VII.3 we presented benchmarks for the Coulomb potential on a lattice that had a dimension of 25.6 Bohr radii. As we found in that section, we were able to determine the ground and first excited states to 0.2% and 2.4%. Improving the accuracy of the first excited state would require going to a lattice with dimensions larger than .
While this is possible with the parallelized code, there is a more efficient way to find excited states by applying symmetry constraints to the initial wavefunction. For example, to find the 1p state of the Coulomb problem we can initialize the wavefunction as as discussed in Sec. IV. In this case we explicitly project out the ground state wavefunction since it is symmetric about the axis. Applying this method on a lattice with lattice spacing and imaginary time step we find the first excited state energy to be which is accurate to 0.06%. At the same time we can extract the next excited state which is antisymmetric about the axis ( state) finding in this case, , corresponding to an accuracy of 0.4%.
The application of symmetry constraints can also allow one to pick out states with different orientations in a 3d potential that breaks spherical symmetry. In Ref Dumitru:2009ni () this technique was used to accurately determine the different heavy quarkonium pwave states corresponding to angular momentum and . Therefore, the ability to constrain the symmetry of the initial FDTD wavefunction is a powerful technique.
vii.7 Application of the Multiresolution Technique
In this section we present benchmarks for the application of the multiresolution technique to the Coulomb potential problem. The current version of the code supports this feature by allowing users the option of saving the wavefunction at the end of the run. The saved wavefunctions can then be read in and used as the initial condition for a subsequent run. The saved wavefunctions can have a different resolution than the resolution of the new run and the code automatically adjusts by sampling/spreading out the wavefunction appropriately.
By using this technique we can accelerate the determination of the high accuracy energy eigenvalues and wavefunctions. In Sec. VII.3 we found that using 32 computational nodes and a random initial wavefunction a run took approximately 1.3 hours. Scaling naively to a lattice, while keeping the lattice volume fixed, would take approximately 42 hours. Using the multiresolution technique and bootstrapping from up to a high resolution ground state and energy eigenvalue can be computed in approximately 45 minutes using the same 32 computational nodes. At the final resolution of and a lattice size of 25.6 Bohr radii the run gives which is accurate to 0.07%. Therefore, the multiresolution technique provides a performance increase of a factor of 50 compared to using random initial wavefunctions for all runs.
Viii Conclusions and Outlook
In this paper we have described a parallel FDTD algorithm for solving the 3d Schrödinger equation. We have shown that for large 3d lattices the method gives a compute time that scales as . This final scaling coefficient and associated error were obtained by averaging the three different scaling coefficients extracted for the Coulomb, harmonic oscillator, and dodecahedron potentials. The crucial optimization that allowed us to achieve nearly ideal scaling was the use of nonblocking sends/receives of the boundary data so that update of each node’s sublattice can proceed while communication of the boundary information is taking place, providing for an “insideout” update algorithm.
Additionally we introduced two novel techniques that can be used in conjunction with the FDTD method. First, we discussed the possibility of imposing a symmetry constraint on the initial wavefunction used for the FDTD evolution. The imposed symmetry constraint allows us to easily construct states that are orthogonal to the ground state and/or some other excited states. Using this technique we can select states that have a certain symmetry, thereby allowing for extremely accurate determination of the particular states we are interested in. Second, we introduced the “multiresolution technique” which simply means that we use the FDTD output wavefunctions from lowerresolution runs as the initial condition for higherresolution runs. Using this method we showed that we can efficiently “bootstrap” our way from small to large lattices, thereby obtaining highaccuracy wavefunctions and eigenvalues in a fraction of the time required when using random initial wavefunctions on large lattices.
The code developed for this paper has been released under an opensource GPL license code (). An obvious next step will be to extend the code to fully 3d partitions, which is in progress. Other areas of improvement include adding support for different types of boundary conditions and complex potentials forth ().
Acknowledgements
We thank V. Antocheviz Dexheimer, A. Dumitru, J. Groff, and J. Milingo for helpful comments. M.S. thanks A. Dumitru, Y. Guo, and A. Mocsy for collaboration in the original project that prompted the development of this code. M.S. also thanks Sharon Stephenson for support of D.Y. during the preparation of this manuscript. D.Y. was supported by the National Science Foundation Award # 0555652.
References
 (1) W.H. Press, et al, Numerical Recipes – The Art of Scientific Computing, Cambridge University Press, 3 edition, Chap. 18 (2007)
 (2) K. S. Yee, IEEE Transacations on Antennas and Propagation 14, 302 (1966).
 (3) A. Taflove, IEEE Trans. Electromagnetic Compatibility 22, 191 (1980).
 (4) A. Taflove and K. R. Umashankar, Electromagnetics 10, 105 (1990).
 (5) A. Taflove, Computational Electrodynamics: The FiniteDifference TimeDomain Method, 3rd edition. Norwood, MA: Artech House (2005).
 (6) V. Varadarajan and R. Mittra, IEEE Microwave and Guided Wave Lett., 4, 144 (1994).
 (7) Z. M. Liu, A.S. Mohan, T.A. Aubrey and W.R Belcher, IEEE Antennas and Propagation Mag. 37, 64 (1995).
 (8) A.D. Tinniswood, P.S. Excell, M. Hargreaves, S. Whittle and D. Spicer, Third International Conference on Computation in Electromagnetics 1, 7 (1996).
 (9) M. Sypniwski, J. Rudnicki and M. CeluchMarcysiak, IEEE Antennas and Propagation International Symposium 1, 252 (2000).
 (10) R. Palaniappan, P. Wahid, G.A. Schiavone and J. Codreanu, IEEE Antennas and Propagation Intern. Symp. 3, 1336 (2000).
 (11) W. Yu, Y. Liu, T. Su and R. Mittra, IEEE Antennas Propag. Mag. 47, 39 (2005).
 (12) W. Yu, M. R. Hashemi, R. Mittra, D. N. de Araujo, N. Pham M. Cases, E. Matoglu, P. Patel, and B. Herrman, IEEE Trans. Adv. Packaging, 30, 335 (2007).
 (13) R. C. Grimm and R. G. Storer, J. Comput. Phys. 7, 134 (1971).
 (14) J. Anderson, J. Chem. Phys. 63, 1499 (1975).
 (15) W. A. Lester Jr., B. L. Hammond, Ann. Rev. Phys. Chem., 41, 283 (1990).
 (16) L. Mitas, Comp. Phys. Commun. 97, 107 (1996).
 (17) M. Foulkes, L. Mitas, R. Needs, G. Rajagopal, Rev. Mod. Phys. 73, 33 (2001).
 (18) N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, N. M. Teller, E. Teller, J. Chem. Phys. 21, 1087 (1953).
 (19) C. Umrigar, Phys. Rev. Lett. 71, 408 (1993).
 (20) D. R. Kent and M. T. Feldmann, QMCBeaver, http://qmcbeaver.sourceforge.net (2004).
 (21) C. Attaccalite, QUMAX, http://qumax.sourceforge.net (2004).
 (22) A. AspuruGuzik et al, J. of Comp. Chem. 26, 856 (2005).
 (23) J.M. Feagin, Quantum Methods with Mathematica, Springer, New York (1994).
 (24) A.K. Roy, N. Gupta and B.M. Deb, Phys. Rev. A 65, 012109 (2001).
 (25) A. Wadehra, A.K. Roy and B.M. Deb, Int. J. Quantum Chem. 91, 597 (2003).
 (26) D.M. Sullivan and D.S. Citrin, J. Appl. Phys. 97, 104305 (2005).
 (27) D.M. Sullivan, J. Appl. Phys. 98, 084311 (2005).
 (28) A.K. Roy, A.J. Thakkar, and B.M. Deb, J. Phys. A: Math. Gen. 38, 2189 (2005).
 (29) I. W. Sudiarta and D. J. W. Geldart, J. Phys. A: Math. Theor. 40, 1885 (2007).
 (30) J. P. Berenger, J. Comput. Phys. 114, 185 (1994).
 (31) W. C. Chew, W. H. Weedon, Microwve Optical Technology Letters 7, 599 (1994).
 (32) S. D. Gedney, IEEE Trans. Antennas Propag. 44, 1630 (1996).
 (33) D. Zhou, W. P. Huang, C. L. Xu, D. G. Fang, and B. Chen, IEEE Photon Technology Letters 13, 454 (2001).
 (34) Y. Rickhard, N. Georgieva, and W. P. Huang, IEEE Microwave Wireless Components Letters 12, 181 (2002).
 (35) O. Ramadan and A. Y. Oztoprak, Microwave Optical Technology Letters 36, 55 (2003).
 (36) O. Ramadan, Parallel Computing 33, 109 (2007).
 (37) O. Ramadan and O. Akaydin, Electr. Eng. 90, 175 (2008).
 (38) T. Fevens and H. Jong, SIAM J. Sci. Comput. 21, 255 (1999).
 (39) I. AlonsoMallo, N. Reguera, Mathematics of Computation 73, 127 (2003).
 (40) H. Han and Z. Xu, Phy. Rev. E 74, 037704 (2004).
 (41) A. R. Mitchell and D. F. Griffiths, The Finite Difference Method in Partial Differential Equations Chichester: John Wiley (1980).
 (42) MPI Forum, MPI: A MessagePassing Interface Standard  Version 2.1, HighPerformance Computing Center Stuttgart, http://www.mpiforum.org/docs/ (2008)

(43)
M. Strickland, Parallelized FDTD Schrödinger Solver,
http://sourceforge.net/projects/quantumfdtd/ (2009).  (44) A. Dumitru, Y. Guo, A. Mocsy and M. Strickland, arXiv:0901.1998 [hepph].
 (45) I. W. Sudiarta and D. J. W. Geldart, Phys. Lett. A, 372(18), 3145 (2008).
 (46) M. Margotta, C. McGahan, D. YagerElorriaga, and M Strickland, forthcoming.