DVR-jl/V_verify.ipynb

88 lines
2.1 KiB
Plaintext

{
"cells": [
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from itertools import chain\n",
"import numpy as np\n",
"import math\n",
"from scipy import sparse\n",
"from scipy.sparse.linalg import eigsh\n",
"\n",
"#sample potential\n",
"def V_gauss(dr_sqr):\n",
" return -4*math.exp(-dr_sqr/4)\n",
"\n",
"n=3 # no of particles\n",
"L=14\n",
"N=6 # no of lattice points\n",
"mu=1/2 # reduced mass\n",
"\n",
"DOF=(n-1)*3 # degrees of freedom after excluding CM\n",
"\n",
"s=np.arange(N**DOF) # matrix index\n",
"\n",
"# k index for each particle and each dimension\n",
"k=np.empty((n-1,3),dtype=np.dtype)\n",
"for dof in range(DOF):\n",
" k[dof//3,dof%3]=(s%N**(DOF-dof))//N**(DOF-1-dof)-N//2\n",
"\n",
"x=k*(L/N) # x coordinate from k index\n",
"\n",
"# adding up all non-local interactions\n",
"V_local=np.zeros(N**DOF)\n",
"\n",
"# 2-body interactions\n",
"V2=np.vectorize(V_gauss)\n",
"dxs=chain((x[i,:] for i in range(n-1)), (x[i,:]-x[j,:] for (i,j) in np.ndindex((n-1,n-1)) if i<j)) # with last particle + with each other\n",
"for dx in dxs:\n",
" dr_sqr=np.sum(dx*dx)\n",
" V_local+=V2(dr_sqr)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"filename = \"temp/V_vals.dat\"\n",
"\n",
"with open(filename, 'br') as f:\n",
" buffer = f.read()\n",
"\n",
"V_julia = np.frombuffer(buffer, dtype=np.float32)\n",
"\n",
"abs_diff = np.abs(V_local-V_julia)\n",
"s = np.sum(abs_diff)\n",
"print(s)"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "base",
"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.13"
},
"orig_nbformat": 4
},
"nbformat": 4,
"nbformat_minor": 2
}