Skip to content

Commit

Permalink
Merge pull request #200 from arnab82/qdpt_bst_pt2
Browse files Browse the repository at this point in the history 
QDPT and SPT-PT2
  • Loading branch information
@nmayhall-vt
nmayhall-vt authored Aug 27, 2024
2 parents c0174ed + 2e42103 commit 7e973f6
Show file tree
Hide file tree
Showing 10 changed files with 1,026 additions and 80 deletions.
22 changes: 8 additions & 14 deletions docs/src/installation_instructions.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,21 +5,15 @@ A Julia package for course-grained electronic structure calculations
### Installation
1. Download

```julia
```
git clone https://github.com/nmayhall-vt/FermiCG.git
cd FermiCG/
cd FermiCG
```


2. Create python virtual environment which will hold the PYSCF executable

```julia
cd src/python
virtualenv -p python3 venv
source venv/bin/activate
pip install -r requirements.txt
cd ../../
julia --project=./
julia> using Pkg; Pkg.build("PyCall")
2. Testing
```
julia --project=./ -tauto
using FermiCG
using Pkg
Pkg.test()
```

327 changes: 327 additions & 0 deletions examples/tetracene_dimer/rhf.ipynb
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,327 @@
{
"cells": [
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Initial guess from minao.\n",
"\n",
"\n",
"******** <class 'pyscf.scf.hf.RHF'> ********\n",
"method = RHF\n",
"initial guess = minao\n",
"damping factor = 0\n",
"level_shift factor = 0\n",
"DIIS = <class 'pyscf.scf.diis.CDIIS'>\n",
"diis_start_cycle = 1\n",
"diis_space = 8\n",
"SCF conv_tol = 1e-09\n",
"SCF conv_tol_grad = None\n",
"SCF max_cycles = 200\n",
"direct_scf = True\n",
"direct_scf_tol = 1e-13\n",
"chkfile to save SCF result = /var/folders/jm/nd85t1m57dv98qh_jtqr2clh0000gn/T/tmpdnsov0fw\n",
"max_memory 4000 MB (current use 0 MB)\n",
"Set gradient conv threshold to 3.16228e-05\n",
"Initial guess from minao.\n",
"init E= -1373.8028556562\n",
" HOMO = -0.155360683388173 LUMO = -0.0462263741041608\n",
"cycle= 1 E= -1360.24317894862 delta_E= 13.6 |g|= 0.719 |ddm|= 11.3\n",
" HOMO = -0.143564650973908 LUMO = 0.138616575609906\n",
"cycle= 2 E= -1360.45923594589 delta_E= -0.216 |g|= 0.214 |ddm|= 1.43\n",
" HOMO = -0.16857409289756 LUMO = 0.12156954691401\n",
"cycle= 3 E= -1360.47370291225 delta_E= -0.0145 |g|= 0.126 |ddm|= 0.515\n",
" HOMO = -0.161685188896432 LUMO = 0.130732383404331\n",
"cycle= 4 E= -1360.4797695373 delta_E= -0.00607 |g|= 0.0139 |ddm|= 0.195\n",
" HOMO = -0.161052562766623 LUMO = 0.132842939686475\n",
"cycle= 5 E= -1360.47997484479 delta_E= -0.000205 |g|= 0.007 |ddm|= 0.0363\n",
" HOMO = -0.161719238589572 LUMO = 0.13333514441787\n",
"cycle= 6 E= -1360.48004918981 delta_E= -7.43e-05 |g|= 0.0028 |ddm|= 0.0275\n",
" HOMO = -0.162116206364872 LUMO = 0.13340292370013\n",
"cycle= 7 E= -1360.48005930686 delta_E= -1.01e-05 |g|= 0.000758 |ddm|= 0.0121\n",
" HOMO = -0.162172259276158 LUMO = 0.133406742943205\n",
"cycle= 8 E= -1360.48005985105 delta_E= -5.44e-07 |g|= 0.000404 |ddm|= 0.00231\n",
" HOMO = -0.162160553391006 LUMO = 0.133425925701453\n",
"cycle= 9 E= -1360.480060043 delta_E= -1.92e-07 |g|= 0.000209 |ddm|= 0.00137\n",
" HOMO = -0.162173349764571 LUMO = 0.133407030042586\n",
"cycle= 10 E= -1360.48006011093 delta_E= -6.79e-08 |g|= 7.63e-05 |ddm|= 0.0011\n",
" HOMO = -0.16217164410856 LUMO = 0.133406153044998\n",
"cycle= 11 E= -1360.48006011748 delta_E= -6.56e-09 |g|= 2.19e-05 |ddm|= 0.000359\n",
" HOMO = -0.162171319204681 LUMO = 0.133405720724885\n",
"cycle= 12 E= -1360.48006011796 delta_E= -4.82e-10 |g|= 6.26e-06 |ddm|= 9.42e-05\n",
" HOMO = -0.162171197951934 LUMO = 0.133405667622305\n",
"Extra cycle E= -1360.480060118 delta_E= -3.27e-11 |g|= 4.12e-06 |ddm|= 1.48e-05\n",
"converged SCF energy = -1360.480060118\n"
]
}
],
"source": [
"from functools import reduce\n",
"import numpy as np\n",
"import scipy\n",
"\n",
"import pyscf\n",
"from pyscf import fci\n",
"from pyscf import gto, scf, ao2mo, lo, tdscf, cc\n",
"\n",
"\n",
"def tda_denisty_matrix(td, state_id):\n",
" '''\n",
" Taking the TDA amplitudes as the CIS coefficients, calculate the density\n",
" matrix (in AO basis) of the excited states\n",
" '''\n",
" cis_t1 = td.xy[state_id][0]\n",
" dm_oo =-np.einsum('ia,ka->ik', cis_t1.conj(), cis_t1)\n",
" dm_vv = np.einsum('ia,ic->ac', cis_t1, cis_t1.conj())\n",
"\n",
" # The ground state density matrix in mo_basis\n",
" mf = td._scf\n",
" dm = np.diag(mf.mo_occ)\n",
"\n",
" # Add CIS contribution\n",
" nocc = cis_t1.shape[0]\n",
" # Note that dm_oo and dm_vv correspond to spin-up contribution. \"*2\" to\n",
" # include the spin-down contribution\n",
" dm[:nocc,:nocc] += dm_oo * 2\n",
" dm[nocc:,nocc:] += dm_vv * 2\n",
"\n",
" # Transform density matrix to AO basis\n",
" mo = mf.mo_coeff\n",
" dm = np.einsum('pi,ij,qj->pq', mo, dm, mo.conj())\n",
" return dm\n",
"\n",
"mol = gto.Mole()\n",
"mol.atom = '''\n",
"C 3.86480 -1.02720 1.27180\n",
"C 3.60170 -0.34020 0.07040\n",
"C 2.93000 -1.06360 2.30730\n",
"C 4.54410 -0.28950 -0.97720\n",
"C 2.32170 0.33090 -0.09430\n",
"C 3.20030 -1.73210 3.54060\n",
"C 1.65210 -0.39240 2.14450\n",
"C 4.27330 0.38740 -2.16700\n",
"C 2.05960 1.01870 -1.29550\n",
"C 1.37970 0.28170 0.95340\n",
"C 2.27260 -1.73270 4.55160\n",
"C 0.71620 -0.41900 3.22490\n",
"C 2.99590 1.05910 -2.32920\n",
"C 5.20900 0.41800 -3.24760\n",
"C 1.01870 -1.06880 4.39470\n",
"C 2.72630 1.73410 -3.55840\n",
"C 4.90670 1.07300 -4.41500\n",
"C 3.65410 1.73950 -4.56900\n",
"C -3.86650 1.01870 -1.29550\n",
"C -2.93020 1.05910 -2.32920\n",
"C -3.60440 0.33090 -0.09430\n",
"C -3.19980 1.73410 -3.55840\n",
"C -1.65280 0.38740 -2.16700\n",
"C -2.32430 -0.34020 0.07040\n",
"C -4.54630 0.28170 0.95340\n",
"C -2.27200 1.73950 -4.56900\n",
"C -0.71710 0.41800 -3.24760\n",
"C -1.38200 -0.28950 -0.97720\n",
"C -2.06130 -1.02720 1.27180\n",
"C -4.27400 -0.39240 2.14450\n",
"C -1.01930 1.07300 -4.41500\n",
"C -2.99610 -1.06360 2.30730\n",
"C -5.20980 -0.41900 3.22490\n",
"C -2.72580 -1.73210 3.54060\n",
"C -4.90730 -1.06880 4.39470\n",
"C -3.65350 -1.73270 4.55160\n",
"H 4.82300 -1.53290 1.39770\n",
"H 5.49910 -0.80290 -0.85660\n",
"H 4.15900 -2.23700 3.66390\n",
"H 1.10180 1.52560 -1.42170\n",
"H 0.42460 0.79440 0.83100\n",
"H 2.50000 -2.24040 5.48840\n",
"H -0.23700 0.09640 3.10140\n",
"H 6.16210 -0.09790 -3.12730\n",
"H 0.29870 -1.07700 5.21470\n",
"H 1.76850 2.24120 -3.67870\n",
"H 5.62580 1.08320 -5.23570\n",
"H 3.42730 2.25190 -5.50340\n",
"H -4.82430 1.52560 -1.42170\n",
"H -4.15760 2.24120 -3.67870\n",
"H -5.50150 0.79440 0.83100\n",
"H -2.49880 2.25190 -5.50340\n",
"H 0.23610 -0.09790 -3.12730\n",
"H -0.42700 -0.80290 -0.85660\n",
"H -1.10300 -1.53290 1.39770\n",
"H -0.30030 1.08320 -5.23570\n",
"H -6.16310 0.09640 3.10140\n",
"H -1.76710 -2.23700 3.66390\n",
"H -5.62740 -1.07700 5.21470\n",
"H -3.42610 -2.24040 5.48840\n",
"'''\n",
"\n",
"\n",
"mol.basis = 'sto-3g'\n",
"mol.spin = 0\n",
"mol.build()\n",
"\n",
"#mf = scf.ROHF(mol).x2c()\n",
"mf = scf.RHF(mol)\n",
"mf.verbose = 4\n",
"mf.get_init_guess(mol, key='minao')\n",
"mf.conv_tol = 1e-9\n",
"#mf.level_shift = .1\n",
"#mf.diis_start_cycle = 4\n",
"#mf.diis_space = 10\n",
"mf.run(max_cycle=200)\n",
"\n",
"\n",
"n_triplets = 2\n",
"n_singlets = 2\n",
"\n",
"avg_rdm1 = mf.make_rdm1()\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
"\n",
"******** <class 'pyscf.tdscf.rhf.TDA'> for <class 'pyscf.scf.hf.RHF'> ********\n",
"nstates = 2 singlet\n",
"deg_eia_thresh = 1.000e-03\n",
"wfnsym = None\n",
"conv_tol = 1e-09\n",
"eigh lindep = 1e-12\n",
"eigh level_shift = 0\n",
"eigh max_space = 50\n",
"eigh max_cycle = 100\n",
"chkfile = /var/folders/jm/nd85t1m57dv98qh_jtqr2clh0000gn/T/tmpdnsov0fw\n",
"max_memory 4000 MB (current use 0 MB)\n",
"\n",
"\n",
"Excited State energies (eV)\n",
"[4.37757462 4.50864348]\n",
"\n",
"** Singlet excitation energies and oscillator strengths **\n",
"Excited State 1: 4.37757 eV 283.23 nm f=0.7064\n",
" 117 -> 123 -0.10590\n",
" 118 -> 124 -0.10581\n",
" 119 -> 122 -0.45930\n",
" 120 -> 121 -0.49153\n",
"Excited State 2: 4.50864 eV 274.99 nm f=0.0000\n",
" 119 -> 121 0.47837\n",
" 120 -> 122 0.46727\n",
"\n",
"** Transition electric dipole moments (AU) **\n",
"state X Y Z Dip. S. Osc.\n",
" 1 2.3095 -1.1112 0.1355 6.5870 0.7064\n",
" 2 -0.0015 0.0006 0.0007 0.0000 0.0000\n",
"\n",
"** Transition velocity dipole moments (imaginary part, AU) **\n",
"state X Y Z Dip. S. Osc.\n",
" 1 -0.1016 0.0712 -0.0291 0.0162 0.0673\n",
" 2 0.0001 0.0000 0.0001 0.0000 0.0000\n",
"\n",
"** Transition magnetic dipole moments (imaginary part, AU) **\n",
"state X Y Z\n",
" 1 -0.0003 -0.0000 0.0022\n",
" 2 -0.0411 -0.1064 0.4355\n"
]
}
],
"source": [
"# compute singlets\n",
"mytd = tdscf.TDA(mf)\n",
"mytd.singlet = True \n",
"mytd = mytd.run(nstates=n_singlets)\n",
"mytd.analyze()\n",
"for i in range(mytd.nroots):\n",
" avg_rdm1 += tda_denisty_matrix(mytd, i)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# compute triplets \n",
"mytd = tdscf.TDA(mf)\n",
"mytd.singlet = False \n",
"mytd = mytd.run(nstates=n_triplets)\n",
"mytd.analyze()\n",
"for i in range(mytd.nroots):\n",
" avg_rdm1 += tda_denisty_matrix(mytd, i)\n",
"\n",
"# normalize\n",
"avg_rdm1 = avg_rdm1 / (n_singlets + n_triplets + 1)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"S = mf.get_ovlp()\n",
"F = mf.get_fock()\n",
"np.save(\"fock_mat18\", F)\n",
"np.save(\"overlap_mat18\", S)\n",
"np.save(\"density_mat18\", mf.make_rdm1())\n",
"np.save(\"rhf_mo_coeffs18\", mf.mo_coeff)\n",
"np.save(\"cis_sa_density_mat18\", avg_rdm1)"
]
}
],
"metadata": {
"kernelspec": {
"display_name": ".venv",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.9.6"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
17 changes: 2 additions & 15 deletions examples/tetracene_dimer/spt_direct.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -39,24 +39,9 @@ FermiCG.add_double_excitons!(ci_vector,FermiCG.FockConfig(init_fspace),nroots)
fspace_0 = FermiCG.FockConfig(init_fspace)
# FermiCG.add_1electron_transfers!(ci_vector, fspace_0, 1)
FermiCG.add_spin_flip_states!(ci_vector, fspace_0,1)
# # Spin-flip states
# ## ba
# tmp_fspace = FermiCG.replace(fspace_0, (1,2), ([4,2],[2,4]))
# FermiCG.add_fockconfig!(ci_vector, tmp_fspace)
# ci_vector[tmp_fspace][FermiCG.TuckerConfig((1:1,1:1))]=FermiCG.Tucker(tuple([zeros(Float64, 1, 1) for _ in 1:nroots]...))
# ## ab
# tmp_fspace = FermiCG.replace(fspace_0, (1,2), ([2,4],[4,2]))
# FermiCG.add_fockconfig!(ci_vector, tmp_fspace)
# ci_vector[tmp_fspace][FermiCG.TuckerConfig((1:1,1:1))]=FermiCG.Tucker(tuple([zeros(Float64, 1, 1) for _ in 1:nroots]...))
display(ci_vector.data)
FermiCG.eye!(ci_vector)
display(ci_vector)
# σ = FermiCG.build_compressed_1st_order_state(ci_vector, cluster_ops, clustered_ham,
# nbody=4,
# thresh=1e-3)
# σ = FermiCG.compress(σ, thresh=1e-5)
# v2 = BSTstate(σ,R=10)
# FermiCG.eye!(v2)
e_ci, v2 = FermiCG.ci_solve(ci_vector, cluster_ops, clustered_ham);
e_var, v_var = FermiCG.block_sparse_tucker(v2, cluster_ops, clustered_ham,
max_iter = 200,
Expand All @@ -70,6 +55,8 @@ e_var, v_var = FermiCG.block_sparse_tucker(v2, cluster_ops, clustered_ham,
resolve_ss = false,
tol_tucker = 1e-4,
solver = "davidson")
@time ept2 = FermiCG.compute_pt2_energy(v_var, cluster_ops, clustered_ham, thresh_foi=1e-6,prescreen = true,compress_twice = true)
@time ept2 = FermiCG.compute_pt2_energy2(v_var, cluster_ops, clustered_ham, thresh_foi=1e-6,prescreen = true,compress_twice = true)
e_var, v_var = FermiCG.block_sparse_tucker(v_var, cluster_ops, clustered_ham,
max_iter = 200,
nbody = 4,
Expand Down
Loading

0 comments on commit 7e973f6

Please sign in to comment.