Doc updates

doc/dchf-mp2.rst

@@ -1,11 +1,11 @@
-*****************************************
-Divide-conquer MP2 by Kobayashi and Nakai
-*****************************************
+*********************************************************
+Divide-conquer perturbation theory by Kobayashi and Nakai
+*********************************************************
 
-Errors in MP2 energies
-======================
+Errors in HF, MP2 and CCSD energies
+===================================
 
-Hartree-Fock (HF) and second-order correction energies (MP2) obtained with DC approach compared to conventional calculations.
+Hartree-Fock (HF) and second-order correction energies (MP2, CCSD) obtained with DC approach compared to conventional calculations.
 
 .. plot:: plots/30-dcmp2-errors_cached.py
 

doc/dchf.rst

@@ -14,7 +14,7 @@ Complexity
 
 The study of how the calculation time scales with the model increase.
 The expected problem complexity in conventional HF is :math:`O(N^4)` while this implementation of DC-HF should scale
-as :math:`O(N^3)` where :math:`N` is the number of atoms in the system.
+as :math:`O(N^2)` where :math:`N` is the number of atoms in the system
 
 .. plot:: plots/21-dchf-complexity_cached.py
 

doc/plots/20-dchf-convergence_cached.py

@@ -3,11 +3,11 @@
 This file was generated automatically.
 """
 from matplotlib import pyplot
-pyplot.semilogy([1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30], [0.0065377003936255562, 0.00016567433431846634, 1.2548197071526879e-05, 9.5038630787414036e-07, 7.1981097304885111e-08, 5.4517605985182627e-09, 4.1291026153800203e-10, 3.1273983402968497e-11, 2.3694379791550091e-12, 1.80855330711438e-13, 1.4543921622589551e-14, 1.7763568394002505e-15, 4.4408920985006262e-16, 7.7715611723760958e-16, 1.1102230246251565e-15, 5.5511151231257827e-16, 6.6613381477509392e-16, 1.1102230246251565e-15, 8.8817841970012523e-16, 6.6613381477509392e-16, 6.6613381477509392e-16, 5.5511151231257827e-16, 1.1102230246251565e-15, 1.1102230246251565e-15, 7.7715611723760958e-16, 8.8817841970012523e-16, 1.1102230246251565e-15, 8.8817841970012523e-16, 6.6613381477509392e-16, 6.6613381477509392e-16], marker='o', label='He chain N=6')
-pyplot.semilogy([1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30], [0.23257339261199045, 0.0064021962035953062, 0.00081903340752387775, 0.00012150868104005008, 1.810213727953397e-05, 2.7807622665354526e-06, 5.8009534542602381e-07, 1.7441402277995621e-07, 5.0205269858105694e-08, 1.3897898637482697e-08, 3.7944926090105469e-09, 1.033047730691905e-09, 2.8175109467731829e-10, 7.7108985863105772e-11, 2.117966912962288e-11, 5.837857974810845e-12, 1.6130985436291212e-12, 4.4639292262615982e-13, 1.247890679678676e-13, 3.5305092183079978e-14, 9.6589403142388619e-15, 3.2474023470285829e-15, 1.0824674490095276e-15, 7.7715611723760958e-16, 1.7486012637846216e-15, 1.3877787807814457e-15, 1.8041124150158794e-15, 7.4940054162198066e-16, 1.2212453270876722e-15, 1.4710455076283324e-15], marker='o', label='H chain N=12 d=6.0')
-pyplot.semilogy([1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30], [0.24272841408793216, 0.011248459301416103, 0.0022875615362965562, 0.00059175210458026206, 0.00016405600492022665, 4.7515044286849939e-05, 1.4517357779769702e-05, 4.3982574390033502e-06, 1.3302581183816109e-06, 4.0271322654383823e-07, 1.2218926009888165e-07, 3.7174423483943286e-08, 1.1342235567424197e-08, 3.4704526197248953e-09, 1.0687561668998313e-09, 3.4305230878217685e-10, 1.0994818944176643e-10, 3.5193875591588153e-11, 1.1253248333176202e-11, 3.5944580645264068e-12, 1.1490253193358058e-12, 3.6698422078984549e-13, 1.1562972801471005e-13, 3.6221026178395732e-14, 1.3711254354120683e-14, 5.4956039718945249e-15, 1.8318679906315083e-15, 2.2204460492503131e-15, 1.5820678100908481e-15, 2.5257573810222311e-15], marker='o', label='H dim chain N=12')
-pyplot.semilogy([1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30], [0.23690318541669264, 0.021185548452057539, 0.0072978632031916856, 0.0026396281061107679, 0.00093216502875580805, 0.0003253849704034327, 0.00011293776323225216, 3.9105704068766178e-05, 1.3525229845305457e-05, 4.6755634414410649e-06, 1.6159824698285322e-06, 5.5850561175331548e-07, 1.9304058162215865e-07, 6.6730812281305674e-08, 2.3071563981340304e-08, 7.9782972164585408e-09, 2.7595025778315296e-09, 9.5464418922025507e-10, 3.3032565745216846e-10, 1.1432075425199528e-10, 3.9572978574597961e-11, 1.3699041900849807e-11, 4.7437210409784214e-12, 1.6459819618397376e-12, 5.6831622741171373e-13, 1.9872124779052314e-13, 6.6787721186845062e-14, 2.3605249699354403e-14, 8.5834117591332415e-15, 3.093011957666647e-15], marker='o', label='+ buffer')
-pyplot.semilogy([1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30], [0.2370892714906648, 0.021628294223540314, 0.0073753621664341389, 0.0026520532129185007, 0.00092460286735897637, 0.00031799794212189567, 0.00010876753084177078, 3.7104518783729994e-05, 1.2641373540937632e-05, 4.3035603160676214e-06, 1.4643475530597486e-06, 4.9807886691975722e-07, 1.6936519083542156e-07, 5.7576429928450634e-08, 1.956939576445299e-08, 6.6502316075842116e-09, 2.2596202348471039e-09, 7.676895612263035e-10, 2.607924804928774e-10, 8.8588626018237804e-11, 3.0091598483827919e-11, 1.0221018476030963e-11, 3.4708512106074174e-12, 1.1800031438080616e-12, 4.0155899438953924e-13, 1.3577680646470469e-13, 4.8057044454985487e-14, 1.5661083541118614e-14, 5.1382509358433026e-15, 2.6645352591003757e-15], marker='o', label='+ larger domain size')
+pyplot.semilogy([1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30], [0.0057099591153211282, 0.00043273689330802778, 3.2776507786658193e-05, 2.482461246233747e-06, 1.8801861201289682e-07, 1.4240299384127297e-08, 1.0785425885728728e-09, 8.1687767661264843e-11, 6.1874949608409224e-12, 4.6918025020659115e-13, 3.6304292905242619e-14, 2.886579864025407e-15, 1.5543122344752192e-15, 6.6613381477509392e-16, 8.8817841970012523e-16, 5.5511151231257827e-16, 7.7715611723760958e-16, 4.4408920985006262e-16, 4.4408920985006262e-16, 6.6613381477509392e-16, 1.1102230246251565e-15, 4.4408920985006262e-16, 8.8817841970012523e-16, 1.3322676295501878e-15, 6.6613381477509392e-16, 9.9920072216264089e-16, 6.6613381477509392e-16, 4.4408920985006262e-16, 8.8817841970012523e-16, 5.5511151231257827e-16], marker='o', label='He chain N=6')
+pyplot.semilogy([1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30], [0.23315574401041875, 0.0054718676826626834, 0.000817068258997692, 0.00012125591082242163, 1.8069102303302012e-05, 2.9965321341474116e-06, 6.3933470637045353e-07, 1.8958003933366996e-07, 5.4476695904304151e-08, 1.5098265720991932e-08, 4.1312619714073406e-09, 1.1274210731659906e-09, 3.0817617902023642e-10, 8.4502876918080005e-11, 2.3249902003641409e-11, 6.4162009039137047e-12, 1.7757462167367066e-12, 4.9260595602618196e-13, 1.3702927681435995e-13, 3.8496983378877303e-14, 1.0380585280245214e-14, 2.4147350785597155e-15, 1.4432899320127035e-15, 7.2164496600635175e-16, 8.6042284408449632e-16, 1.0547118733938987e-15, 1.4710455076283324e-15, 1.3322676295501878e-15, 7.4940054162198066e-16, 8.6042284408449632e-16], marker='o', label='H chain N=12 d=6.0')
+pyplot.semilogy([1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30], [0.2444390577770294, 0.011219809219615656, 0.0027152449123842848, 0.00071754503757975874, 0.00020133409996225815, 6.0686633627876629e-05, 1.8383412336647753e-05, 5.5335516817578423e-06, 1.6640562366943001e-06, 5.010539270988712e-07, 1.5122903121445042e-07, 4.5772065804250417e-08, 1.3894643990930433e-08, 4.2303346103889083e-09, 1.2916762381909308e-09, 3.9549444275088774e-10, 1.2389445025462464e-10, 3.9709902033280287e-11, 1.271105443123588e-11, 4.0641656706696949e-12, 1.2984058272991206e-12, 4.1436298836572405e-13, 1.3233858453531866e-13, 4.2577052994374753e-14, 1.3600232051658168e-14, 6.0229599085914742e-15, 2.5257573810222311e-15, 1.3045120539345589e-15, 1.6653345369377348e-15, 1.4988010832439613e-15], marker='o', label='H dim chain N=12')
+pyplot.semilogy([1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30], [0.23730106384430072, 0.019232158047645789, 0.0069718986250189158, 0.0025213243799438317, 0.00088841774242279931, 0.00030989044185442821, 0.00010749415839633042, 3.720098002045133e-05, 1.2860113957128411e-05, 4.4435563134326284e-06, 1.5350965525198768e-06, 5.3031316482946012e-07, 1.8321516776984015e-07, 6.3306546487645887e-08, 2.1878107291095628e-08, 7.5623052403495095e-09, 2.6144914219705662e-09, 9.0408685746168338e-10, 3.1269623348983977e-10, 1.0817488918513085e-10, 3.7428775356840305e-11, 1.2952881822680951e-11, 4.4844995772397311e-12, 1.5519634188887466e-12, 5.3780244146928169e-13, 1.8668747103767203e-13, 6.3640065439685145e-14, 2.1732615707037439e-14, 7.5564554613549717e-15, 3.4416913763379853e-15], marker='o', label='+ buffer')
+pyplot.semilogy([1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30], [0.23745459996264384, 0.019651223410304633, 0.007038587717119963, 0.0025290833609339786, 0.00088126031571042138, 0.00030285996744594881, 0.00010349762395550632, 3.5271658253746563e-05, 1.2004244052440938e-05, 4.0821416095217222e-06, 1.3874143690306706e-06, 4.7135251342658724e-07, 1.6008214048096603e-07, 5.4352671297747124e-08, 1.8450088359878158e-08, 6.261679154351274e-09, 2.1247686908532071e-09, 7.2089711977696425e-10, 2.4456089237068213e-10, 8.295878220487829e-11, 2.8138779500919142e-11, 9.5448891845872907e-12, 3.2372472758002147e-12, 1.0984720077988896e-12, 3.7305575295576432e-13, 1.2576398256136656e-13, 4.25735835474228e-14, 1.4988010832439613e-14, 5.3117232834409833e-15, 4.163336342344337e-15], marker='o', label='+ larger domain size')
 pyplot.xlabel('Step')
 pyplot.ylabel('Maximal absolute change in amplitudes')
 pyplot.grid(axis='y')

doc/plots/21-dchf-complexity.py

@@ -34,14 +34,24 @@ def setup_full(model):
     return hf
 
 
+def get_scaling(x, y):
+    while True:
+        try:
+            return curve_fit(lambda x, a, p: a * (x ** p), x, y, p0=(5, y[-1]/x[-1]**5))[0]
+        except RuntimeError:
+            x, y = x[1:], y[1:]
+            if len(x)<2:
+                return float("nan")
+
+
 for label, constructor in (
         ("He chains", helium_chain),
         ("H chains", hydrogen_dimer_chain),
 ):
-    for driver, label2, marker in (
-        (scf.RHF, "pyscf HF", "x"),
-        (setup, "DCHF", "o"),
-        (setup_full, "custom HF", "+")
+    for driver, label2, marker, ls in (
+        (scf.RHF, "pyscf HF", "x", "--"),
+        (setup, "DCHF", "o", None),
+        (setup_full, "custom HF", "+", "--")
     ):
         dt = None
 
@@ -63,17 +73,14 @@ def setup_full(model):
                 times.append(dt)
                 print "  time {:.3f}s".format(times[-1])
 
-        f = lambda x, a, p: a*(x**p)
-        try:
-            popt = curve_fit(f, chain_size[:len(times)], times, p0=(5, times[-1]/chain_size[len(times)-1]**5))[0]
-        except RuntimeError:
-            popt = [float("nan"), float("nan")]
+        popt = get_scaling(chain_size[:len(times)][-3:], times[-3:])
 
         pyplot.loglog(
             chain_size[:len(times)],
             times,
             label=label+" ("+label2+") s={:.1f}".format(popt[1]),
             marker=marker,
+            ls=ls,
         )
 
 pyplot.xlabel("Model size")

doc/plots/21-dchf-complexity_cached.py

@@ -3,12 +3,12 @@
 This file was generated automatically.
 """
 from matplotlib import pyplot
-pyplot.loglog([8, 12, 16, 24, 32, 48, 64], [0.22306489944458008, 0.2650279998779297, 0.5159330368041992, 0.9762320518493652, 1.7665889263153076, 7.157909870147705, 23.577895879745483], marker='x', label='He chains (pyscf HF) s=4.0')
-pyplot.loglog([8, 12, 16, 24, 32, 48], [0.2416379451751709, 0.4748649597167969, 0.947746992111206, 2.9169440269470215, 6.480128049850464, 20.54787302017212], marker='o', label='He chains (DCHF) s=2.8')
-pyplot.loglog([8, 12, 16, 24, 32], [0.18092012405395508, 0.38437581062316895, 1.0169260501861572, 4.117337942123413, 12.579202890396118], marker='+', label='He chains (custom HF) s=3.8')
-pyplot.loglog([8, 12, 16, 24, 32, 48], [0.3520781993865967, 0.41600680351257324, 0.6842679977416992, 1.3300771713256836, 2.995187997817993, 11.026674032211304], marker='x', label='H chains (pyscf HF) s=3.1')
-pyplot.loglog([8, 12, 16, 24, 32], [0.5784149169921875, 1.248039960861206, 2.0974950790405273, 5.964249849319458, 12.80048394203186], marker='o', label='H chains (DCHF) s=2.6')
-pyplot.loglog([8, 12, 16, 24, 32], [0.3278329372406006, 0.9977419376373291, 2.206886053085327, 8.124358892440796, 22.903008937835693], marker='+', label='H chains (custom HF) s=3.5')
+pyplot.loglog([8, 12, 16, 24, 32, 48, 64], [0.31583499908447266, 0.3838682174682617, 0.3919641971588135, 0.9799468517303467, 1.8762309551239014, 8.902273893356323, 25.437674045562744], marker='x', ls='--', label='He chains (pyscf HF) s=3.7')
+pyplot.loglog([8, 12, 16, 24, 32, 48, 64], [0.2991189956665039, 0.6342098712921143, 0.710867166519165, 1.5551469326019287, 3.8979990482330322, 9.137217044830322, 18.05000901222229], marker='o', ls=None, label='He chains (DCHF) s=2.3')
+pyplot.loglog([8, 12, 16, 24, 32, 48], [0.19790410995483398, 0.42900705337524414, 0.8198988437652588, 2.966587781906128, 8.431477069854736, 41.83495283126831], marker='+', ls='--', label='He chains (custom HF) s=3.9')
+pyplot.loglog([8, 12, 16, 24, 32, 48, 64], [0.18721699714660645, 0.5095930099487305, 0.6442949771881104, 1.2542729377746582, 2.255833864212036, 9.29486608505249, 29.024532079696655], marker='x', ls='--', label='H chains (pyscf HF) s=3.9')
+pyplot.loglog([8, 12, 16, 24, 32, 48], [0.44788503646850586, 0.8136281967163086, 1.282900094985962, 3.3723721504211426, 6.198496103286743, 13.524383068084717], marker='o', ls=None, label='H chains (DCHF) s=2.0')
+pyplot.loglog([8, 12, 16, 24, 32], [0.27031588554382324, 0.6370410919189453, 1.5198769569396973, 5.466639995574951, 15.520428895950317], marker='+', ls='--', label='H chains (custom HF) s=3.6')
 pyplot.xlabel('Model size')
 pyplot.ylabel('Run time (s)')
 pyplot.grid()

doc/plots/22-dchf-errors.py

@@ -1,7 +1,7 @@
 from pyscf import scf
 
 import dchf
-from test_common import helium_chain, hydrogen_dimer_chain, hydrogen_distant_dimer_chain
+from test_common import helium_chain, hydrogen_dimer_chain, hydrogen_distant_dimer_chain, atomic_chain
 from test_dchf import assign_domains
 
 import fake
@@ -13,34 +13,46 @@ def calculate(model, domain_size, buffer_size=0):
     assign_domains(hf, domain_size, buffer_size)
 
     hf.kernel(tolerance=1e-12)
-    return hf.hf_energy
+    return hf.hf_energy, hf.dm
 
 models = [
     (helium_chain(6), dict(domain_size=1), "He chain N=6"),
     (hydrogen_distant_dimer_chain(12), dict(domain_size=2), "H chain N=12 d=6.0"),
     (hydrogen_dimer_chain(12), dict(domain_size=2), "H dim chain N=12"),
-    (hydrogen_dimer_chain(12), dict(domain_size=2, buffer_size=2), "+ buffer"),
-    (hydrogen_dimer_chain(12), dict(domain_size=4, buffer_size=2), "+ larger domain size"),
+    (hydrogen_dimer_chain(12), dict(domain_size=4, buffer_size=2), "... 3(0) config"),
+    (hydrogen_dimer_chain(12), dict(domain_size=2, buffer_size=1), "... 2(1) config"),
+    (hydrogen_dimer_chain(12), dict(domain_size=2, buffer_size=2), "... 2(2) config"),
+    (hydrogen_dimer_chain(12), dict(domain_size=4, buffer_size=2), "... 4(2) config"),
 ]
 
 
 names = []
 e = []
+dm = []
 
 for model, options, label in models:
     print label
     mf = scf.RHF(model)
     mf.conv_tol = 1e-12
     mf.kernel()
-    e_tot = calculate(model, **options)
+    e_tot, hf_dm = calculate(model, **options)
     e.append(abs(mf.e_tot - model.energy_nuc() - e_tot)/model.natm)
+    dm.append(abs(hf_dm - mf.make_rdm1()).max())
     names.append(label)
 
 x = list(range(len(names)))
-pyplot.figure(figsize=(8, 4.8))
-pyplot.semilogy(x, e, marker="_", ls="None", markersize=30, mew=5)
+pyplot.figure(figsize=(16, 4.8))
 
-pyplot.xticks(x, names)
+pyplot.subplot(121)
+pyplot.semilogy(x, e, marker="_", ls="None", markersize=30, mew=5)
+pyplot.xticks(x, names, rotation=15)
 pyplot.ylabel("Energy difference per atom (Ry)")
 pyplot.grid(axis='y')
+
+pyplot.subplot(122)
+pyplot.semilogy(x, dm, marker="_", ls="None", markersize=30, mew=5)
+pyplot.xticks(x, names, rotation=15)
+pyplot.ylabel("Maximal density matrix deviation")
+pyplot.grid(axis='y')
+
 pyplot.show()

doc/plots/22-dchf-errors_cached.py

@@ -3,10 +3,16 @@
 This file was generated automatically.
 """
 from matplotlib import pyplot
-pyplot.figure(figsize=(8, 4.8))
-pyplot.semilogy([0, 1, 2, 3, 4], [4.2040445199139263e-14, 2.1260164739800302e-08, 0.0016855996096232968, 3.4306672232660894e-05, 1.5154790137176613e-05], marker='_', mew=5, markersize=30, ls='None')
-pyplot.xticks([0, 1, 2, 3, 4], ['He chain N=6', 'H chain N=12 d=6.0', 'H dim chain N=12', '+ buffer', '+ larger domain size'])
+pyplot.figure(figsize=(16, 4.8))
+pyplot.subplot(121)
+pyplot.semilogy([0, 1, 2, 3, 4, 5, 6], [1.0835776720341528e-13, 2.126016444374083e-08, 0.001685599609623889, 1.4123599610874985e-05, 0.011493718776240947, 3.4306672231032564e-05, 1.4123599610578927e-05], marker='_', mew=5, markersize=30, ls='None')
+pyplot.xticks([0, 1, 2, 3, 4, 5, 6], ['He chain N=6', 'H chain N=12 d=6.0', 'H dim chain N=12', '... 3(0) config', '... 2(1) config', '... 2(2) config', '... 4(2) config'], rotation=15)
 pyplot.ylabel('Energy difference per atom (Ry)')
 pyplot.grid(axis='y')
+pyplot.subplot(122)
+pyplot.semilogy([0, 1, 2, 3, 4, 5, 6], [1.7481654457363049e-08, 0.00027446799566078865, 0.023748974905231514, 0.0027834788210318526, 0.23184446461122366, 0.005309024815691385, 0.0027834788210312073], marker='_', mew=5, markersize=30, ls='None')
+pyplot.xticks([0, 1, 2, 3, 4, 5, 6], ['He chain N=6', 'H chain N=12 d=6.0', 'H dim chain N=12', '... 3(0) config', '... 2(1) config', '... 2(2) config', '... 4(2) config'], rotation=15)
+pyplot.ylabel('Maximal density matrix deviation')
+pyplot.grid(axis='y')
 pyplot.show()
 

doc/plots/23-dchf-errors-domain.py

@@ -24,7 +24,7 @@
 pyplot.subplot(121)
 for domain_size in [2, 4, 6, 8, 12, 24]:
     print "Size={:d}".format(domain_size)
-    x = list(sorted(set([0, 2, (domain_size//4)*2, domain_size, 2*domain_size])))
+    x = list(sorted(set([0, 2, (domain_size//4)*2, domain_size])))
     y = []
     for buffer_size in x:
         print "  buff={:d}".format(buffer_size)