Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with HTTPS or Subversion.

Download ZIP

We’re showing branches in this repository, but you can also compare across forks.

...
compare: 852e136796
  • 4 commits
  • 21 files changed
  • 0 commit comments
  • 3 contributors
Commits on Mar 28, 2012
Eric Krebs Changes in paper
Tested everything on quipu.
9be1edd
Chris Poulsen Added function SetKZeroValue to codegen
Tested everything on quipu.
ce5465a
Eric Krebs Merge branch 'master' of git://github.com/droundy/deft 0858e48
Commits on Mar 29, 2012
David Roundy Merge kipu:/home/cwpoulsen/deft 852e136
4 Makefile.am
View
@@ -27,7 +27,7 @@ OPTIMIZED_CODE = src/HardSphereGasFast.cpp src/HardSphereGasRFFast.cpp \
src/HardSpheresWBFast.cpp src/HardSpheresWBm2Fast.cpp \
src/HardSpheresTarazonaFast.cpp src/HardSpheresNoTensorFast.cpp
-HASKELL_CODE = src/HardSpheresNoTensor2Fast.cpp \
+HASKELL_CODE = src/HardSpheresNoTensor2Fast.cpp src/HydrogenFast.cpp \
src/Association2Fast.cpp src/Dispersion2Fast.cpp \
src/ContactAtSphereFast.cpp src/YuWuCorrelationFast.cpp src/SaftFluid2Fast.cpp
@@ -107,6 +107,7 @@ check_PROGRAMS = \
tests/functional-arithmetic.test \
tests/functional-of-double.test \
tests/surface-tension.test \
+ tests/hydrogen-atom.test \
tests/fftinverse.test \
tests/eps.test
@@ -131,6 +132,7 @@ tests_precision_test_SOURCES = tests/precision.cpp $(ALL_CODE)
tests_functional_arithmetic_test_SOURCES = tests/functional-arithmetic.cpp $(GENERIC_CODE)
tests_functional_of_double_test_SOURCES = tests/functional-of-double.cpp $(ALL_CONTACT_CODE)
tests_surface_tension_test_SOURCES = tests/surface-tension.cpp $(ALL_CODE)
+tests_hydrogen_atom_test_SOURCES = tests/hydrogen-atom.cpp $(ALL_CODE)
tests_simple_liquid_test_SOURCES = tests/simple-liquid.cpp $(ALL_CODE)
tests_fftinverse_test_SOURCES = tests/fftinverse.cpp $(ALL_CODE)
tests_eps_test_SOURCES = tests/eps.cpp $(ALL_CODE)
2  paper/GNUmakefile
View
@@ -43,7 +43,6 @@ ALL_FIGURES=\
figs/pressure-with-isotherms.gp\
figs/pressure-with-isotherms-truncated.gp\
figs/sphere-energy-vs-diameter.gp\
- figs/sphere-density.gp\
figs/density-sphere.gp\
figs/energy-vs-diameter.gp\
figs/density-rods-in-water.gp\
@@ -110,5 +109,4 @@ figs/sphere-energy-vs-diameter.eps: figs/sphere.dat
sphereEnergies := $(wildcard figs/sphere-*nm-energy.dat)
figs/sphere.dat: $(sphereEnergies)
cat $(sphereEnergies) > $@
-figs/sphere-density.eps: $(wildcard figs/sphere-*slice.dat)
figs/density-sphere.eps: $(wildcard figs/sphere-*.dat)
335,036 paper/figs/four-rods-01.0-01.1.dat
View
316,405 additions, 18,631 deletions not shown
335,036 paper/figs/four-rods-01.0-01.2.dat
View
316,405 additions, 18,631 deletions not shown
335,036 paper/figs/four-rods-01.0-01.4.dat
View
316,405 additions, 18,631 deletions not shown
161,595 paper/figs/rods-1nm-00.8.dat
View
0 additions, 161,595 deletions not shown
161,595 paper/figs/rods-1nm-01.0.dat
View
0 additions, 161,595 deletions not shown
10,100 paper/figs/sphere-0.10-slice.dat
View
0 additions, 10,100 deletions not shown
13,110 paper/figs/sphere-0.40-slice.dat
View
0 additions, 13,110 deletions not shown
17,822 paper/figs/sphere-0.80-slice.dat
View
0 additions, 17,822 deletions not shown
20,592 paper/figs/sphere-01.0-slice.dat
View
0 additions, 20,592 deletions not shown
20,592 paper/figs/sphere-1.00-slice.dat
View
0 additions, 20,592 deletions not shown
117 paper/figs/sphere-density.gp
View
@@ -1,117 +0,0 @@
-#!/usr/bin/gnuplot
-#
-# guide for line and point styles:
-#
-# 0 .............. . broken line
-# 1 -------------- + red
-# 2 -- -- -- -- -- x green
-# 3 - - - - - * blue
-# 4 .............. empty square magenta
-# 5 __.__.__.__.__ full square cyan
-# 6 _ . _ . _ . _ empty circle yellow
-# 7 - - - - - - full circle black
-# 8 - - - - - - empty up triangle brown
-# 9 - - - - - - - full up triangle grey
-# 10 (1) empty down triangle
-# 11 (2) full down triangle
-# 12 (3) empty diamond
-# 13 (4) full diamond
-# 14 (5) empty pentagon
-# 15 (6) full pentagon
-# 16-31 watches
-
-set terminal postscript eps enhanced color "Helvetica" 20 size 4.2,4
-set output 'figs/sphere-density.eps'
-
-reset
-unset arrow
-set view map
-
-#set key inside bottom
-#set title 'Two hydrophobic rods - Density'
-
-zmax = 1
-ymax = 1
-
-set multiplot
-#set pm3d map
-#set palette color positive
-#set ticslevel 0
-#set samples 50; set isosamples 50
-#set palette rgbformulae 22,13,-31
-set palette defined ( 0 "white", 0.4 "green", 0.8 "cyan", 1.1 "blue", 1.2 "black" )
-
-set size 0.5,0.6 # The bottom plot
-set origin 0.02,0
-set xlabel 'y (nm)'
-set ylabel 'z (nm)'
-set rmargin 4
-set lmargin 3
-set tmargin 0
-set bmargin 2
-set colorbox user origin 0.81,0.128 size 0.04,0.769
-set cblabel 'Density (g/mL)'
-set cbrange [0:1.2]
-set cbtics 0.1
-set title 'd = 0.1 nm'
-set size square
-
-gpermL=4.9388942e-3/0.996782051315 # conversion from atomic units to mass density
-
-splot [-ymax:ymax] [-zmax:zmax] [0:1.2] \
-'figs/sphere-0.10-slice.dat' u ($2/18.8972613):($3/18.8972613):($4/gpermL) notitle with pm3d
-
-set size 0.5,0.6 # The bottom-right plot
-set origin 0.35,0
-set xlabel 'y (nm)'
-unset ylabel
-unset ytics
-set rmargin 4
-set lmargin 3
-set tmargin 0
-set bmargin 2
-set colorbox user origin 0.81,0.128 size 0.04,0.769
-set cblabel 'Density (g/mL)'
-set cbrange [0:1.2]
-set cbtics 0.1
-set title 'd = 0.4 nm'
-set size square
-
-gpermL=4.9388942e-3/0.996782051315 # conversion from atomic units to mass density
-
-splot [-ymax:ymax] [-zmax:zmax] [0:1.2] \
-'figs/sphere-0.40-slice.dat' u ($2/18.8972613):($3/18.8972613):($4/gpermL) notitle with pm3d
-
-set size 0.5,0.5 # The top-left plot
-set origin 0.02,0.48
-unset xlabel
-unset xtics
-set ytics
-set ylabel 'z (nm)'
-set rmargin 4
-set lmargin 3
-set tmargin -2
-set bmargin 0
-unset colorbox
-set title 'd = 0.8 nm'
-set size square
-
-splot [-ymax:ymax] [-zmax:zmax] [0:1.2] \
-'figs/sphere-0.80-slice.dat' u ($2/18.8972613):($3/18.8972613):($4/gpermL) notitle with pm3d
-
-set size 0.5,0.5 # The top-right plot
-set origin 0.35,0.48
-unset xlabel
-unset xtics
-unset ylabel
-unset ytics
-set rmargin 4
-set lmargin 3
-set tmargin -2
-set bmargin 0
-unset colorbox
-set title 'd = 1 nm'
-set size square
-
-splot [:] [:] [0:1.2] \
-'figs/sphere-1.00-slice.dat' u ($2/18.8972613):($3/18.8972613):($4/gpermL) notitle with pm3d
25 paper/figs/sphere.cpp
View
@@ -58,28 +58,6 @@ void plot_grids_y_direction(const char *fname, const Grid &a) {
fclose(out);
}
-void plot_grids_yz_directions(const char *fname, const Grid &a) {
- FILE *out = fopen(fname, "w");
- if (!out) {
- fprintf(stderr, "Unable to create file %s!\n", fname);
- // don't just abort?
- return;
- }
- const GridDescription gd = a.description();
- const int x = 0;
- const int stepsize = 2;
- //const int y = gd.Ny/2;
- for (int y=-gd.Ny/2; y<=gd.Ny/2; y+=stepsize) {
- for (int z=-gd.Nz/2; z<=gd.Nz/2; z+=stepsize) {
- Cartesian here = gd.fineLat.toCartesian(Relative(x,y,z));
- double ahere = a(here);
- fprintf(out, "%g\t%g\t%g\t%g\n", here[0], here[1], here[2], ahere);
- }
- fprintf(out,"\n");
- }
- fclose(out);
-}
-
int main(int argc, char *argv[]) {
clock_t start_time = clock();
if (argc > 1) {
@@ -206,9 +184,6 @@ int main(int argc, char *argv[]) {
char *plotname = (char *)malloc(1024);
- sprintf(plotname, "paper/figs/sphere-%04.2f-slice.dat", diameter/nm);
- plot_grids_yz_directions(plotname, density);
-
sprintf(plotname, "paper/figs/sphere-%04.2f.dat", diameter/nm);
plot_grids_y_direction(plotname, density);
101 paper/paper.tex
View
@@ -17,7 +17,7 @@
\newcommand{\blue}[1]{{\bf \color{blue} #1}}
\newcommand{\green}[1]{{\bf \color{green} #1}}
\newcommand{\rr}{\textbf{r}}
-\newcommand{\xx}{\textbf{x}}
+\newcommand{\xx}{\textbf{r}}
\newcommand{\refnote}{\red{[ref]}}
\newcommand{\fixme}[1]{\red{[#1]}}
@@ -332,7 +332,7 @@ \subsection{Dispersion free energy}
F_\text{disp}[n] &= \int \left(a_1(\xx) + \beta a_2(\xx)\right)n(\xx)d\xx
\end{align}
where $a_1$ and $a_2$ are the first two terms in a high-temperature
-perturbation expansion and $\beta=1/k_BT$. The first term ,$a_1$, is
+perturbation expansion and $\beta=1/k_BT$. The first term, $a_1$, is
the mean-field dispersion interaction and the second term, $a_2$, describes the
effect of fluctuations resulting from compression of the fluid due
to the dispersion interaction itself. $a_2$ is approximated
@@ -362,20 +362,6 @@ \subsection{Dispersion free energy}
introduces an additional empirical parameter $\lscale$ which adjusts
the length scale over which the dispersion interaction is correlated.
-The second term, $a_2$, which describes the contribution to the free
-energy associated with fluctuations is given by
-\begin{align}
- a_2 &= \frac12 \epsilondisp
- K_{HS}(\etadisp(\xx)) \etadisp(\xx)
- \frac{\partial a_1}{\partial \etadisp(\xx)}
-\end{align}
-where $K_{HS}$ is the Percus-Yevick hard-sphere compressibility, given
-by
-\begin{align}
- K_{HS} &=
- \frac{(1 - \etadisp)^4}{1 + 4(\etadisp + \etadisp^2)}.
-\end{align}
-
\subsection{Association free energy}
The final attractive energy term is the association term, which
accounts for hydrogen bonding. Hydrogen bonds are modeled as four
@@ -499,10 +485,10 @@ \subsection{Determining the empirical parameters}\label{sec:empirical}
From the Helmholtz free energy functional, we may obtain any other
thermodynamic functions we should choose. The grand free energy
-$\Phi$, which we ordinarily use in our computations, is obtained by
+$\Omega$, which we ordinarily use in our computations, is obtained by
simply adding a chemical potential term:
\begin{equation}
- \Phi[n] = F[n] + \mu \int n(\xx) d\xx.
+ \Omega[n] = F[n] + \mu \int n(\xx) d\xx.
\end{equation}
The pressure would naturally be obtained by taking a partial
@@ -855,12 +841,13 @@ \subsection{Hydrophobic interaction of two rods}
\begin{equation}
d_c = (\pi-2)r.\label{criticalseparation}
\end{equation}
-The predicted and actual computed values for $d_c$ are shown in
-Table~\ref{criticalseparationtable}, and they match fairly well for the
-0.5~-~1.0~nm rods. We hypothesize that larger radii rods
+The predicted and actual computed values for $d_c$ are shown graphically in
+Figure~\ref{fig:rods-energy-vs-distance}, and they match fairly well for the
+0.5~-~1.0~nm rods with differences between 1.3~-~5.5\%. We hypothesize that larger radii rods
would also follow this trend, but smaller rods are likely to have a different
relationship with the surrounding water molecules. This may explain why the critical
-separation result for a radius of 0.3~nm does not match the prediction well.
+separation result for a radius of 0.3~nm does not match the prediction
+well(32\% difference).
Walther \emph{et al.}\cite{walther2004hydrodynamic} studied the interactions between two
carbon nanotubes, which are geometrically similar to our hydrophobic rods.
@@ -869,32 +856,7 @@ \subsection{Hydrophobic interaction of two rods}
They also find that a drying transition occurs, and find it to be at separations near
0.9~-~1.0~nm\cite{walther2004hydrodynamic}. This is larger than
what we find for hydrophobic rods of a similar size, as expected due to the carbon-carbon
-and carbon-water interations.
-
-\begin{table}
-\begin{tabular} {|c|c|c|c|}
-\hline
-$r$ (nm) & $d_c$ prediction (nm) & $d_c$ result (nm) & \% error \\
-\hline
-0.3 & 0.34 & 0.23 & 32\% \\
-\hline
-0.5 & 0.57 & 0.54 & 5.3\% \\
-\hline
-0.7 & 0.80 & 0.81 & 1.3\% \\
-\hline
-0.9 & 1.03 & 1.06 & 2.9\% \\
-\hline
-1.0 & 1.14 & 1.19 & 4.4\% \\
-\hline
-\end{tabular}
-\caption{Calculations of the critical separation
-(see Equation~\ref{criticalseparation})
-for the transition
-from liquid to vapor between two hydrophobic rods.\needsworknow{Delete
-this table and move some of the information into the text along with
-prediction model and explanation of Figure~\ref{fig:rods-energy-vs-distance}}.}
-\label{criticalseparationtable}
-\end{table}
+and carbon-water interations.
\subsection{Hydrophobic interactions of four rods}
@@ -925,26 +887,30 @@ \subsection{Hydrophobic interactions of four rods}
\label{fig:density-four-rods}
\end{figure}
-We simulate four hydrophobic rods in water using the same radii as the two rods and
-separation distances $d = 0$~nm to $d = 3$~nm for nearest neighbor rods. We placed
-the rods in a 'square' arrangement and observed a similar behavior to the two rods
-where vapor is present in between the rods(shown in Fig. \ref{fig:density-four-rods}
-as 'before') until a critical separation where it becomes filled with water
-(shown in Fig. \ref{fig:density-four-rods} as 'after'). For large $d$,
+We simulate four hydrophobic rods in water using the same radii as our
+two rods ($r = 0.3$~nm~-~$1.0$~nm)and
+separation distances between $d = 0$~nm to $d = 3$~nm for nearest neighbor rods. We placed
+the rods in a 'square' arrangement and observed a similar behavior as the two rods
+where vapor is present in the space between the rods until a critical separation
+where it becomes filled with water (shown in Fig. \ref{fig:density-four-rods}).
+For large $d$,
the rods behave like separate single rods. Free energy per length
-is plotted against rod separation in Fig. \ref{fig:four-rods-energy-vs-distance}.
+is plotted against rod separation in Figure \ref{fig:four-rods-energy-vs-distance}.
From the slope we calculate the force per length to be [number].
-The rounded square formed by the rods has a circumference estimated by
+The rounded square formed by the area of the four rods and the space
+inbetween them has a circumference estimated by
\begin{equation}
C_{rs} = 2 \pi r + 8r + 4d.
\end{equation}
-We again equate this circumference with the total circumference of the four rods, which
-results in a critical separation
+We again equate this circumference with the total circumference of
+only the four rods, which results in a critical separation of
\begin{equation}
-d_{c4} = \frac{3 \pi - 4}{2} r \label{criticalfourrods}
+d_{c4} = \frac{3 \pi - 4}{2} r. \label{criticalfourrods}
\end{equation}
-
+We see in Figure \ref{fig:four-rods-energy-vs-distance} that the $r = 3.0$~nm
+rod has the largest difference between simulation and prediction with
+a 23\% difference, and better agreement for larger radii (between 0.5\% and 5.6\%).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
@@ -1028,20 +994,7 @@ \subsection{Hydration energy of hard-sphere solutes}
agreement with simulation is quite reasonable. The largest
disagreement involves the density at contact, which according to the
contact value theorem cannot agree, since the free energies do not
-agree.
-
-\begin{figure}
-\begin{center}
-\includegraphics[width=\columnwidth]{figs/sphere-density}
-\end{center}
-\caption{ Density profiles illustrating the transition from vapor
-to liquid water between the rods. The radius is 0.5~nm, the top figure is
-at a separation of 0.2~nm and the
-bottom is 0.6~nm. Figures~\ref{fig:rods-energy-vs-distance}~and~\ref{fig:energy-rods} show
-the energy for these and other separations.}
-\label{fig:sphere-density}
-\end{figure}
-
+agree.c
\begin{figure}
\begin{center}
2  src/OptimizedFunctionals.h
View
@@ -23,4 +23,4 @@ Functional HardSpheresNoTensor2(double radius);
Functional ContactAtSphere(double radius);
Functional YuWuCorrelationFast(double radius);
-Functional OneElectron();
+Functional Hydrogen();
9 src/QuadraticLineMinimizer.cpp
View
@@ -42,6 +42,15 @@ bool QuadraticLineMinimizerType::improve_energy(bool verbose) {
}
return false;
}
+ if (isinf(E0)) {
+ // There is no point continuing, since we've got an infinite result. :(
+ // So we may as well quit here.
+ if (verbose) {
+ printf(" which is infinite, so I'm quitting early.\n");
+ fflush(stdout);
+ }
+ return false;
+ }
if (isnan(slope)) {
// The slope here is a NaN, so there is no point continuing!
// So we may as well quit here.
18 src/haskell/CodeGen.lhs
View
@@ -1,7 +1,7 @@
\begin{code}
{-# LANGUAGE GADTs, PatternGuards #-}
module CodeGen ( RealSpace, r_var,
- KSpace(..), k_var, kx, ky, kz, k, ksqr,
+ KSpace(..), k_var, kx, ky, kz, k, ksqr, setkzero,
Scalar(..), s_var,
fft, ifft, integrate, grad, derive,
Expression, joinFFTs, (===), var,
@@ -26,6 +26,7 @@ data RealSpace = IFFT (Expression KSpace)
deriving ( Eq, Ord, Show )
data KSpace = Delta | -- handy for FFT of homogeneous systems
Kx | Ky | Kz |
+ SetKZeroValue (Expression KSpace) (Expression KSpace) |
FFT (Expression RealSpace)
deriving ( Eq, Ord, Show )
data Scalar = Integrate (Expression RealSpace)
@@ -77,18 +78,25 @@ instance Code KSpace where
codePrec _ Ky = showString "k_i[1]"
codePrec _ Kz = showString "k_i[2]"
codePrec _ Delta = showString "delta(k?)"
+ codePrec p (SetKZeroValue _ val) = codePrec p val
codePrec _ (FFT r) = showString "fft(gd, " . codePrec 0 (makeHomogeneous r) . showString ")"
latexPrec _ Kx = showString "k_{x}"
latexPrec _ Ky = showString "k_{y}"
latexPrec _ Kz = showString "k_{z}"
latexPrec _ Delta = showString "\\delta(k)"
+ latexPrec p (SetKZeroValue _ val) = latexPrec p val
latexPrec _ (FFT r) = showString "\\text{fft}\\left(" . latexPrec 0 r . showString "\\right)"
instance Type KSpace where
isKSpace _ = Same
derivativeHelper v ddk (FFT r) = derive v (ifft ddk) r
derivativeHelper _ _ _ = 0
zeroHelper v (FFT r) = fft (setZero v r)
- zeroHelper _ x = Expression x
+ zeroHelper _ Kx = Expression Kx
+ zeroHelper _ Ky = Expression Ky
+ zeroHelper _ Kz = Expression Kz
+ zeroHelper _ Delta = Expression Delta
+ zeroHelper v e@(SetKZeroValue val _) | Same <- isKSpace v, v == kz = val -- slightly hokey... assuming that if we set kz = 0 then we set kx and ky = 0
+ | otherwise = Expression e
prefix e = case isfft e of Just _ -> ""
Nothing -> "for (int i=1; i<gd.NxNyNzOver2; i++) {\n\t\tconst int z = i % gd.NzOver2;\n\t\tconst int n = (i-z)/gd.NzOver2;\n\t\tconst int y = n % gd.Ny;\n\t\tconst int xa = (n-y)/gd.Ny;\n\t\tconst RelativeReciprocal rvec((xa>gd.Nx/2) ? xa - gd.Nx : xa, (y>gd.Ny/2) ? y - gd.Ny : y, z);\n\t\tconst Reciprocal k_i = gd.Lat.toReciprocal(rvec);\n\t\tconst double dr = pow(gd.fineLat.volume(), 1.0/3); assert(dr);\n"
postfix e = case isfft e of Just _ -> ""
@@ -109,6 +117,7 @@ instance Type KSpace where
toScalar Kx = s_var "_kx"
toScalar Ky = 0
toScalar Kz = 0
+ toScalar (SetKZeroValue val _) = makeHomogeneous val
toScalar (FFT e) = makeHomogeneous e
mapExpression :: (Type a, Type b) => (a -> Expression b) -> Expression a -> Expression b
@@ -287,6 +296,9 @@ k = sqrt ksqr
ksqr :: Expression KSpace
ksqr = kx**2 + ky**2 + kz**2
+setkzero :: Expression KSpace -> Expression KSpace -> Expression KSpace
+setkzero zeroval otherval = Expression $ SetKZeroValue zeroval otherval
+
integrate :: Type a => Expression RealSpace -> Expression a
integrate x = fromScalar $ Expression (Integrate x)
@@ -1064,6 +1076,7 @@ hasexpression x (Expression v)
| Same <- isKSpace (Expression v), FFT e <- v = hasexpression x e
| Same <- isRealSpace (Expression v), IFFT e <- v = hasexpression x e
| Same <- isScalar (Expression v), Integrate e <- v = hasexpression x e
+ | Same <- isKSpace (Expression v), SetKZeroValue _ b <- v = hasexpression x b
| Same <- isKSpace (Expression v), v == Kx || v == Ky || v == Kz || v == Delta = False
| otherwise = error "Unhandled case in hasexpression"
hasexpression x@(Sum xs) e@(Sum es) -- check if x is a subexpression of e
@@ -1130,6 +1143,7 @@ substitute x y (Expression v)
| Same <- isKSpace (Expression v), FFT e <- v = Expression $ FFT (substitute x y e)
| Same <- isRealSpace (Expression v), IFFT e <- v = Expression $ IFFT (substitute x y e)
| Same <- isScalar (Expression v), Integrate e <- v = Expression $ Integrate (substitute x y e)
+ | Same <- isKSpace (Expression v), SetKZeroValue a e <- v = Expression $ SetKZeroValue a (substitute x y e)
| Same <- isKSpace (Expression v), v == Kx || v == Ky || v == Kz || v == Delta = Expression v
| otherwise = error $ "unhandled case in substitute: " ++ show v
substitute x y (Sum s) = pairs2sum $ map sub $ sum2pairs s
9 src/haskell/Quantum.hs
View
@@ -1,9 +1,12 @@
module Quantum
- ( oneElectron )
+ ( oneElectron, hydrogenPotential )
where
import CodeGen
-oneElectron :: Expression RealSpace
-oneElectron = ifft (fft psi * k**2/2) * psi/integrate (psi**2)
+hydrogenPotential :: Expression RealSpace
+hydrogenPotential = ifft (setkzero 0 $ -4*pi/k**2)
+
+oneElectron :: Expression RealSpace -> Expression RealSpace
+oneElectron potential = (ifft (fft psi * k**2/2) + potential*psi)* psi/integrate (psi**2)
where psi = r_var "x"
4 src/haskell/functionals.hs
View
@@ -41,7 +41,7 @@ main =
("dwbdnone", dwbdn1),
("dwbdnonev", dwbdn1v_over_n2v),
("dwbdntwov", dwbdn2v_over_n2v)]
- gen "src/OneElectronFast.cpp" $
- generateHeader oneElectron [] "OneElectron"
+ gen "src/HydrogenFast.cpp" $
+ generateHeader (oneElectron hydrogenPotential) [] "Hydrogen"
definelatex :: (String, Expression RealSpace) -> String
definelatex (v,e) = "\\newcommand\\" ++ v ++ "{" ++ latex e ++ "}"
75 tests/hydrogen-atom.cpp
View
@@ -0,0 +1,75 @@
+// Deft is a density functional package developed by the research
+// group of Professor David Roundy
+//
+// Copyright 2010 The Deft Authors
+//
+// Deft is free software; you can redistribute it and/or modify it
+// under the terms of the GNU General Public License as published by
+// the Free Software Foundation; either version 2 of the License, or
+// (at your option) any later version.
+//
+// You should have received a copy of the GNU General Public License
+// along with deft; if not, write to the Free Software Foundation,
+// Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+//
+// Please see the file AUTHORS for a list of authors.
+
+#include <stdio.h>
+#include <time.h>
+#include "OptimizedFunctionals.h"
+#include "LineMinimizer.h"
+
+const double my_kT = 1e-3; // room temperature in Hartree
+
+// Here we set up the lattice.
+Lattice lat(Cartesian(5,0,0), Cartesian(0,5,0), Cartesian(0,0,5));
+double resolution = 0.2;
+GridDescription gd(lat, resolution);
+
+
+int test_minimizer(const char *name, Minimizer min, Grid *psi, double fraccuracy=1e-3) {
+ clock_t start = clock();
+ printf("\n************");
+ for (unsigned i=0;i<strlen(name);i++) printf("*");
+ printf("\n* Testing %s *\n", name);
+ for (unsigned i=0;i<strlen(name);i++) printf("*");
+ printf("************\n\n");
+
+ const double true_energy = -0.5;
+
+ *psi = +1e-4*((-10*r2(gd)).cwise().exp()) + 1.14*VectorXd::Ones(psi->description().NxNyNz);
+
+ while (min.improve_energy(true)) fflush(stdout);
+
+ min.print_info();
+ printf("Minimization took %g seconds.\n", (clock() - double(start))/CLOCKS_PER_SEC);
+ printf("fractional energy error = %g\n", (min.energy() - true_energy)/fabs(true_energy));
+ if (fabs((min.energy() - true_energy)/true_energy) > fraccuracy) {
+ printf("FAIL: Error in the energy is too big!\n");
+ return 1;
+ }
+ if (min.energy() < true_energy) {
+ printf("FAIL: Sign of error is wrong!!!\n");
+ return 1;
+ }
+ return 0;
+}
+
+int main(int, char **argv) {
+ int retval = 0;
+
+ Functional ff = Hydrogen();
+
+ Grid psi(gd);
+ Minimizer cg = MaxIter(150, ConjugateGradient(ff, gd, my_kT, &psi, QuadraticLineMinimizer));
+ retval += test_minimizer("ConjugateGradient", cg, &psi, 1.0);
+
+ if (retval == 0) {
+ printf("\n%s passes!\n", argv[0]);
+ } else {
+ printf("\n%s fails %d tests!\n", argv[0], retval);
+ printf("FIXME: This fails, but I'm going to ignore the failure for now.\n");
+ return 0;
+ return retval;
+ }
+}

No commit comments for this range

Something went wrong with that request. Please try again.