From f4d390972de4da50340283ea0e3dc5961f632f3f Mon Sep 17 00:00:00 2001 From: Panadestein Date: Mon, 26 Apr 2021 21:50:11 +0200 Subject: [PATCH] Solutions in C up to 4 --- .gitignore | 1 + hydrogen.h | 14 ++++++++ pp | Bin 0 -> 17136 bytes qmc_gaussian.c | 7 ++++ qmc_stats.h | 34 ++++++++++++++++++ tt.c | 11 ++++++ vmc_metropolis.c | 87 +++++++++++++++++++++++++++++++++++++++++++++++ 7 files changed, 154 insertions(+) create mode 100755 pp create mode 100644 qmc_gaussian.c create mode 100644 tt.c create mode 100644 vmc_metropolis.c diff --git a/.gitignore b/.gitignore index c912f89..e7ebbb4 100644 --- a/.gitignore +++ b/.gitignore @@ -5,3 +5,4 @@ get_energ get_var qmc_uni qmc_met +vmc diff --git a/hydrogen.h b/hydrogen.h index 905b482..c061238 100644 --- a/hydrogen.h +++ b/hydrogen.h @@ -46,3 +46,17 @@ double kinetic(double a, double *r, const int l) { double e_loc(double a, double *r, const int l) { return kinetic(a, r, l) + potential(r, l); } + +void drift(double a, double *r, double *d, const int l) { + double rnorm, fact; + + for (int i = 0; i < l; ++i) { + rnorm += r[i]*r[i]; + } + rnorm = sqrt(rnorm); + + fact = -a / rnorm; + for (int i = 0; i < l; ++i) { + d[i] = r[i] * fact; + } +} diff --git a/pp b/pp new file mode 100755 index 0000000000000000000000000000000000000000..30dfb73be233038c1fa369d67cbf9ca5159d5d3e GIT binary patch literal 17136 zcmeHPaco=1c|VGh%``DdC-$N^P0Xe-DrXQaJB};Ib1lkx_ERUvt}V4TV&|C_DT#}c?K1ShZfm>}-tatrBD z1#MeY2Bbr@i)Hw`THGs^fNzwTRo-C{lv?R>LnZYo+zU#2&6KG^KWM>%DSJqi^lGI& zssg4WXi1V@6J-_erYmo;L0VJkQJ&MiYR#_|#Z>^-~R9vPMWj$g^ zZ&c}xDm|uxy`Z&TkM)(c59 zqTNM%Em$z+`aS_YipzgDaX0T(_0{T!vnq~ERd;UIqV54g5?Ed`}Jh zKn?uQYv7O9z&`>!gv*_d0;uNyp&Ixm;2~V@MB}WQo#7hz9|OP7ZG?D?mn033O~UBz z-DSk$>G)tGW5(0HyS5LfM&iBEBg1iF7=y{wh>u>B5jKPR7PYH zBV@ndvY$yuM`B_$Yi5L*NRmyM>wnznA9~yvh$e=`XgV=s4gl{z7BvPEBhlf+lO#qa zRMp8SOz!OI?%Zx{2yX~Kw26<4jA7d-romGp`6wvKQ(BQoETd8XBjgFQ^c$o;mGivG^VMP1Swm`GWIUvB zUk7O|%x}zJ-e%!+50Sr|3J!cJgP=tRt{x4t^o9e!%^`o&fz$oU=?e-c+Y~)c6Aql` z3@UrZf%Du)_*n-opJN0fJ?p@&`N6`^IdF~MY0iQB zZM2@LN1z^odIahbs7Ih4f&U8;_(kB}xAol5eEL|^kJkyIpPn`w%5!?|g73U^zx>c% zz@_qp#pmlI$QAft~r?Ip7LrnwSb=chcn?v`>Pu=*sTA<-&lB6YEp|p4>{sBKy02dT#6Ul-2Xs%-i&_tuKNp z=HDxqi!lWM#is9r_Z)(*T|ep1e;S39>4V-rjJ9??f8ugJdr8mt&E`&A@?_WOW7_3X zASgg*OM&H*>bYs(K;VsO;^&uRj9t_Us9Ribv!sH)F|D9KqmQq6fnYMBd*;YWDX@Qr zE#Zc71u+)(4D`4I+|LraUnP6 z#iOsO zQ?H2|<+7gtSAFfPdcLG@zGU96pZXbX_4M4C`)>vyetjO2r(P=~udl_J;QZ-V%-|Y{ z=88VuQ_yohGzDaOsXkNm>N5o|oG$z=RarFODVOsXi-$liYA+FVy`}g?a0nwtP;l86 z>}ph(7smQt()(XQ2pZ>d+Dqt($t~JTCuBGZp3~ET(|5|zHcMva-^L$AA4;=OLQLln z{*me8E59b=!DX*w8$~@#9poI0Q~LM^8;b8hio0vEgdDCIBg^}qMdA1!4^<>x3pQdT zAJFECpMinuiU+i3p@}|n*DI|yqwG*kt9`ar8UPD-(3HG`l;xZ@+W_lZw8_Akw{ert ze;;ANgffTuw)husQ=~7rLUf@bMBj!Ch3Fy_tPowO5~3eNY9T~!3E7ws8WTrjE>sQC zei%TR%QPw4A*xmwY=PBeO{*JhtG_Ws7Z9RJ+&`PO3xP9lAxhMRV||kqf%wb+v;y%1 zJh`kuOuD?Ku`&m#g+RC^WFXq$auzO6Rt>~%7(kiJv?CBb`KuG~uNc!Gw5A5LsrdR$ z_&grJ@o41W!c2l`q|G^h+{FAbt`&YjPjVgi*&ky3T-2^m?_a-vOuK?A$GTo9$U6fg zWD3u(DOwd2zY6O)?TRON;!49)uTGJIKGt^=4y=%K6C6lD% z-T7A|eckzAMS3Ioclz}4)=9wode**6YnI}+x8VBByIynk`k$!hT~EH$lmB^F{@)_y z6+hK;=RNx7f6l&6>zPLnMGi&|MIMV7=LZJDfz_WQU*!6LRys6=lFfyck1VeLsYjq5 zfqDe$5vWJt|78T|_ZUjup=dHRgdeTM(dck!&G3NwaiF=itM}I*x~l^}Z(Aut#|Q_a`_9O&tR#44)j@2T084+mdm4{>r3Tw0d)Kq<#Hod1^Q zr@y(=AMEt| z->0cqlo=o+}&KR=qg3|`#9QvP{MXareD2RPEqD8 zy~#!5y@IYX{0-TxD(7#+{?QT_{B2mX;#VrlM5&2U^?jJNyA*`(J+0(uze$Pfzt>{z_Tjy>Y#-t(Rq;Yz@*T<^HnEbo{ZjnBsvK|7WjVJ) zng1Un_IKPai*2gEI~3ic=wU@i6+Nlwgra8^J*VijqAv6A%lbMjk)7MOKOAc5JCYqS zv!Tu52g7Zx>pz;6?E0rRw1wL?uC?IzwMbuedhGA}o#OO8DSw|{iF*Zqe_n|l631O#rC+=_MG>ls z`^8*UJRk~H@t}QPqq5%Rg6EA&93RVuoQRe8iaYG@@tvh!yd{~aD!ys}Q> zb_1vOy6ry-JhWI`WNX;@yu^dzg*z-I`FLWQi95>^!aC_UPw(3H3v?F0e)9L7O zBR*oLkBfnHG#NKy*<|uKR2-ZENwY%KG~tk8?AROGr5V~IU33_wQYt2lt^##HJELKfx30_zcu5 zkMS6>Ov)I-u^c+$)4dx~vBZdx&BSB!h>qKm9TDml7(0^5D8uqxPV{Kph^N!3G#&hL z>*6>O_ZOWrs(N^c4h~i7(UGB4(in_pGc@F!2ZJgFML2UjX-1ELnrVv;aSo@7;^|Ql z9!Z(;@Zd-`Jep38#?$6;2XrKxz?r#3OhJ*(?p8B8D5UsMG&3Z^vEw7C!lGu{DmfZY zXA-HA3dTTLIzAjF1(h2eHbq#*ER1Y;Fa?$}nRveloAF~{>7IaeN)DTFd`R5|Lorxm z&az?M4weSzQB^dVfH6xQ{UXA+6_dEVs@&SY4MZyxOci+cS*cmAYnk$XO%+R7&J+DG zF4`lpKCh#hhEV2AynkYh&cQkL_pADu((h@`#LorBx{-0}^SYdAhy@jO>wgM3ehN#m z|Ga)@%Jr-5n7yxKvr%Md56b$y&S!d(f{R2xQkL7V9xEaNjGD~)ydPl7`&@{~q|7lZ z=-edjAz7dI4NOBylKp2nrca@q_N2`7{(|XYrOz>D{}n98k)fEfKJP=APAfsKpY^%@ zUs3v41fl*A@J|5;)Qq;WyclZBSUHT!V$Ml>s$HZk`e8Z*B z`y{3;UrV3tJN+*x1*T_NQBk(X^f{M4uWy*1R`Fyrtj9ckt?I0w_h(G`eGkjK<45an z(&zpcAtkA3yJDGe|1tk}$T;=eLskjXUKUi;UH?yjJN0>A$aD`CI+DBoYbc{VmYR0# z{Ue{xcFQx*^iAjn9QwS^9BKuK#Bvt5z%%|X3TW(d{k*^7^Y$EA?nc(*ex&^pHPO~r z@3}_R0(Sys&cymmmq6L6FFKTTn?<=Ew!wPL-v(vExPD&e29>^Bp6YTY+D8$l#O)&1 zf51kp4C_0oJ@y>39m-%x1@0b{;lP32A3Xoi`#bjCp)L;7MjDC#FM`EczH*-Ea$ME) SgAZF$tJqPCm2Fj7N&X9IDaUvK literal 0 HcmV?d00001 diff --git a/qmc_gaussian.c b/qmc_gaussian.c new file mode 100644 index 0000000..1d537df --- /dev/null +++ b/qmc_gaussian.c @@ -0,0 +1,7 @@ +#include +#include + +double gaussian(double *r) { + // Pending + return 0; +} diff --git a/qmc_stats.h b/qmc_stats.h index 3ac432d..0bd615f 100644 --- a/qmc_stats.h +++ b/qmc_stats.h @@ -1,5 +1,7 @@ #include #include +#include +#include void ave_error(double *x, const int n, double obs[2]) { // obs[1] --> mean and obs[2] --> err @@ -25,3 +27,35 @@ void ave_error(double *x, const int n, double obs[2]) { } } + +void random_gauss(double *z, const int n) { + // Box Muller method + const double two_pi = 2.0 * acos(-1.0); + double u[n + 1]; + + srand(time(NULL)); + + for (int i = 0; i < n + 1; ++i) { + u[i] = (double) rand() / RAND_MAX; + } + + + if (!(n & 1)) { + // n is even + for (int i = 0; i < n; i+=2) { + z[i] = sqrt(-2.0 * log(u[i])); + z[i] = z[i] * cos(two_pi * u[i + 1]); + z[i + 1] = z[i] * sin(two_pi * u[i + 1]); + } + } + else { + // n is odd + for (int i = 0; i < n - 1; i+=2) { + z[i] = sqrt(-2.0 * log(u[i])); + z[i] = z[i] * cos(two_pi * u[i + 1]); + z[i + 1] = z[i] * sin(two_pi * u[i + 1]); + } + z[n] = sqrt(-2.0 * log(u[n])); + z[n] = z[n] * cos(two_pi * u[n + 1]); + } +} diff --git a/tt.c b/tt.c new file mode 100644 index 0000000..20f4d12 --- /dev/null +++ b/tt.c @@ -0,0 +1,11 @@ +#include "qmc_stats.h" + +int main() { + double rnd[3]; + random_gauss(rnd, 3); + for (int i = 0; i < 3; ++i) { + printf("val %lf\n", rnd[i]); + } + + return 0; +} diff --git a/vmc_metropolis.c b/vmc_metropolis.c new file mode 100644 index 0000000..2175e34 --- /dev/null +++ b/vmc_metropolis.c @@ -0,0 +1,87 @@ +#include "hydrogen.h" +#include "qmc_stats.h" + +void variational_montecarlo(double a, const int nmax, double dt, + double *energy, double *accept) { + int n_accept; + double r_old[3], r_new[3], psi_old, psi_new; + double d_old[3], d_new[3], d2_old, d2_new; + double rnd[3], aval, fact_a, fact_b, fact_exp, u; + + // Initial position + random_gauss(r_old, 3); + psi_old = psi(a, r_old, 3) * psi(a, r_old, 3); + drift(a, r_old, d_old, 3); + d2_old = 0.0; + for (int i = 0; i < 3; ++i) { + d2_old += d_old[i] * d_old[i]; + } + + *energy = 0.0; + *accept = 0.0; + n_accept = 0; + for (int i = 0; i < nmax; ++i) { + // Compute and accumulate the local energy + *energy += e_loc(a, r_old, 3); + + // Compute new position + random_gauss(rnd, 3); + for (int j = 0; j < 3; ++j) { + r_new[j] = r_old[j] + dt * d_old[j] + rnd[j]; + } + + // New WF and acceptance probability + psi_new = psi(a, r_new, 3) * psi(a, r_new, 3); + + drift(a, r_new, d_new, 3); + d2_new = 0.0; + for (int i = 0; i < 3; ++i) { + d2_new += d_new[i] * d_new[i]; + } + + // Compute the ratio of probabilities q + fact_a = 0.5 * dt * (d2_new - d2_old); + fact_b = 0.0; + for (int i = 0; i < 3; ++i) { + fact_b += (d_new[i] + d_old[i]) * (r_new[i] - r_old[i]); + } + fact_exp = exp(fact_b - fact_a); + aval = fact_exp * psi_new / psi_old; + + u = (double) rand() / RAND_MAX; + + if (u <= aval) { + for (int j = 0; j < 3; ++j) { + r_old[j] = r_new[j]; + d_old[j] = d_new[j]; + } + psi_old = psi_new; + n_accept += 1; + } + } + *energy /= nmax; + *accept = (double) n_accept / nmax; +} + +int main() { + const double a = 1.2; + const double dt = 1.0; + const long nmax = 1e5; + int nruns = 30; + + srand(time(NULL)); + + double ene[nruns], acc[nruns], obs[2]; + + for (int i = 0; i < nruns; ++i) { + variational_montecarlo(a, nmax, dt, &ene[i], &acc[i]); + } + + ave_error(ene, nruns, obs); + printf("E = %.5lf +/- %.5lf\n", obs[0], obs[1]); + + ave_error(acc, nruns, obs); + printf("A = %.5lf +/- %.5lf\n", obs[0], obs[1]); + + return 0; +}