diff --git a/el_nuc_el.irp.f b/el_nuc_el.irp.f index 84cffd4..2f342ee 100644 --- a/el_nuc_el.irp.f +++ b/el_nuc_el.irp.f @@ -3,7 +3,7 @@ BEGIN_PROVIDER [double precision, factor_een] BEGIN_DOC ! Electron-electron nucleus contribution to Jastrow factor END_DOC - integer :: i, j, alpha, p, k, l, lmax + integer :: i, j, alpha, p, k, l, lmax, cindex double precision :: x, y, z, t, c_inv, u, a, b, a2, b2, c, t0 PROVIDE cord_vect @@ -19,6 +19,7 @@ BEGIN_PROVIDER [double precision, factor_een] b2 = b * b c = rescale_een_n(i, alpha) * rescale_een_n(j, alpha) c_inv = 1.0d0 / c + cindex = 0 do p = 2, ncord x = 1.0d0 do k = 0, p - 1 @@ -47,7 +48,8 @@ BEGIN_PROVIDER [double precision, factor_een] ! This can be used in case of a flatten cord_vect ! cidx = 1 + l + (ncord + 1) * k + (ncord + 1) * (ncord + 1) * (p - 1) + & ! (ncord + 1) * (ncord + 1) * ncord * (alpha - 1) - factor_een = factor_een + cord_vect(l, k, p, alpha) * (y + z) * t + cindex = cindex + 1 + factor_een = factor_een + cord_vect(cindex, typenuc_arr(alpha)) * (y + z) * t t = t * c_inv y = y * a2 z = z * b2 diff --git a/elec_coord.txt b/elec_coord.txt index 0d33424..50e632f 100644 --- a/elec_coord.txt +++ b/elec_coord.txt @@ -1,10 +1,10 @@ - 0.47107089936082125 0.36934601008972345 0.13813129163530469 - 0.11734479025488309 0.36584274529092520 0.92672173271006275 - 0.35754797210216693 0.24703841579308339 0.17042939824115788 - 0.31813480208419698 0.81855922566341688 0.93785308019180225 - 0.69675380000569165 0.48537277058577866 2.8246759780116260E-002 - 0.42646034508932040 5.1852406712776689E-002 0.93402091856761427 - 0.85418292218111624 0.51831660694819415 0.26738545106704503 - 0.34322893096420737 0.80132196716714288 0.33887619652274625 - 0.69754618540855695 0.39691148679374866 0.12716214781484358 - 0.57090264858892203 0.70538842967387250 0.31300897404191264 + -0.250655104764153 0.503070975550133 -0.166554344502303 + -0.587812193472177 -0.128751981129274 0.187773606533075 + 1.61335569047166 -0.615556732874863 -1.43165470979934 + -4.901239896295210E-003 -1.120440036458986E-002 1.99761909330422 + 0.766647499681200 -0.293515395797937 3.66454589201239 + -0.127732483187947 -0.138975497694196 -8.669850480215846E-002 + -0.232271834949124 -1.059321673434182E-002 -0.504862241464867 + 1.09360863531826 -2.036103063808752E-003 -2.702796910818986E-002 + -0.108090166832043 0.189161729653261 2.15398313919894 + 0.397978144318712 -0.254277292595981 2.54553335476344 diff --git a/electrons.irp.f b/electrons.irp.f index 0a7bed7..b83b836 100644 --- a/electrons.irp.f +++ b/electrons.irp.f @@ -13,14 +13,15 @@ BEGIN_PROVIDER [ double precision, elec_coord, (nelec, 3) ] ! Electron coordinates END_DOC character(len=*), parameter :: FILE_NAME = "elec_coord.txt" - integer :: fu, rc + integer :: fu, rc, i, j open(action='read', file=FILE_NAME, iostat=rc, newunit=fu) - read(fu, *) elec_coord + do i = 1, nelec + read(fu, *) elec_coord(i, :) + end do close(fu) - END_PROVIDER BEGIN_PROVIDER [ double precision, elec_dist, (nelec, nelec) ] @@ -28,15 +29,16 @@ BEGIN_PROVIDER [ double precision, elec_dist, (nelec, nelec) ] BEGIN_DOC ! e-e distance END_DOC - integer :: i,j - double precision :: x,y,z + integer :: i, j + double precision :: x, y, z + do j = 1, nelec - do i = 1, nelec - x = elec_coord(i, 1) - elec_coord(j, 1) - y = elec_coord(i, 2) - elec_coord(j, 2) - z = elec_coord(i, 3) - elec_coord(j, 3) - elec_dist(i, j) = dsqrt( x*x + y*y + z*z ) - enddo + do i = 1, nelec + x = elec_coord(i, 1) - elec_coord(j, 1) + y = elec_coord(i, 2) - elec_coord(j, 2) + z = elec_coord(i, 3) - elec_coord(j, 3) + elec_dist(i, j) = dsqrt( x*x + y*y + z*z ) + enddo enddo END_PROVIDER @@ -47,18 +49,19 @@ BEGIN_PROVIDER [double precision, factor_ee] END_DOC integer :: i, j, p double precision :: pow_ser, x + factor_ee = 0.0d0 pow_ser = 0.0d0 - do j = 1 , nelec + do j = 1, nelec do i = 1, nelec - x = rescale_ee(i, j) + x = rescale_ee(i, j) * rescale_ee(i, j) do p = 2, nbord x = x * rescale_ee(i, j) - pow_ser = pow_ser + bord_vect(p) * x + pow_ser = pow_ser + bord_vect(p + 1) * x end do factor_ee = factor_ee + bord_vect(1) * rescale_ee(i, j) & - / (1 + bord_vect(2) * rescale_ee(i, j)) + pow_ser + / (1.0d0 + bord_vect(2) * rescale_ee(i, j)) + pow_ser end do end do diff --git a/jast_coeffs.txt b/jast_coeffs.txt index 593d63b..a064ab3 100644 --- a/jast_coeffs.txt +++ b/jast_coeffs.txt @@ -1,190 +1,35 @@ - 4.7107089936082125E-004 - 1.1734479025488309E-004 - 3.5754797210216693E-004 - 3.1813480208419699E-004 - 6.9675380000569162E-004 - 4.7107089936082121E-008 - 1.1734479025488308E-008 - 3.5754797210216693E-008 - 3.1813480208419699E-008 - 6.9675380000569162E-008 - 4.7107089936082132E-006 - 1.1734479025488310E-006 - 3.5754797210216696E-006 - 3.1813480208419699E-006 - 6.9675380000569168E-006 - 4.2646034508932041E-006 - 8.5418292218111638E-006 - 3.4322893096420741E-006 - 6.9754618540855701E-006 - 5.7090264858892204E-006 - 3.6934601008972349E-006 - 3.6584274529092522E-006 - 2.4703841579308343E-006 - 8.1855922566341692E-006 - 4.8537277058577867E-006 - 5.1852406712776692E-007 - 5.1831660694819417E-006 - 8.0132196716714300E-006 - 3.9691148679374872E-006 - 7.0538842967387255E-006 - 1.3813129163530470E-006 - 9.2672173271006290E-006 - 1.7042939824115789E-006 - 9.3785308019180231E-006 - 2.8246759780116260E-007 - 9.3402091856761429E-006 - 2.6738545106704506E-006 - 3.3887619652274627E-006 - 1.2716214781484360E-006 - 3.1300897404191267E-006 - 2.5932160037276120E-006 - 8.9781735169153588E-006 - 4.4307945771984676E-006 - 3.0357326911046815E-006 - 4.7057217460270767E-006 - 9.1888685615387477E-006 - 1.3674140578661999E-006 - 2.0898016957168522E-006 - 8.3133097492563700E-006 - 8.1930838386368430E-006 - 5.0761544134895837E-006 - 9.6316690786857319E-006 - 7.3265410239880617E-006 - 8.0282246083863469E-006 - 1.7840082668643953E-006 - 8.7402085917097983E-007 - 5.0021478344850869E-006 - 6.9080963589020322E-006 - 3.6461563934505928E-006 - 2.8079367396118341E-006 - 6.4561942747179404E-006 - 2.0143276334222562E-006 - 7.2282413910517550E-006 - 4.2803612481228009E-006 - 9.8791262118408863E-006 - 1.0651141027935296E-008 - 2.1266859661958340E-006 - 4.7254708764414140E-006 - 6.7960239791661172E-006 - 2.5361974009558710E-006 - 8.3022917885344452E-006 - 3.9304102414302272E-007 - 4.6287756683245075E-006 - 4.8818348901487870E-006 - 5.4217534412207837E-006 - 9.2515623629987977E-006 - 5.6776909044188313E-008 - 8.5143283852193907E-006 - 8.4916543326313185E-006 - 2.6803664649533946E-006 - 3.8619825277037558E-007 - 4.1645240139124852E-007 - 4.7076790048846678E-006 - 3.7400607062400029E-006 - 6.5903674693304497E-007 - 7.0511416830091020E-006 - 7.1970472254075054E-006 - 8.3779145190230263E-006 - 4.4720525973186892E-006 - 5.3527028821450744E-006 - 3.2469611963900162E-006 - 2.7474685177806738E-006 - 2.8766396048376057E-006 - 3.0832590779442896E-007 - 1.1970794212078929E-006 - 9.8875410686779822E-006 - 3.8833096533945047E-006 - 6.5493531824651745E-006 - 4.8669577823889157E-006 - 3.0536887907279188E-006 - 1.5453212737038092E-006 - 4.2010386366569512E-006 - 3.4608479854625175E-006 - 5.5863755844856395E-006 - 3.1626318255858814E-006 - 8.9535689383160287E-006 - 2.7295097461499133E-006 - 6.5761403856116881E-006 - 9.9777710552222946E-006 - 5.7483205746228453E-006 - 7.2938932476699272E-006 - 1.8975608710106486E-006 - 9.0217837038189077E-006 - 9.2351027820514674E-006 - 2.6491154807706522E-006 - 1.5921008779667090E-006 - 8.6057966599747186E-006 - 3.7446440900687352E-006 - 3.4272640538904732E-006 - 2.0533004904915601E-006 - 9.2668744161565381E-006 - 5.7347639398840740E-006 - 2.1537537756204565E-006 - 4.9515420151455622E-006 - 2.5237493063925278E-006 - 8.7662988099822665E-006 - 6.5819974241574931E-007 - 9.2929324306096999E-006 - 3.4693265321361314E-006 - 4.5904275891192394E-007 - 4.3467474285524594E-006 - 5.6775174141423421E-007 - 4.2394412605675084E-006 - 6.9012994065703607E-008 - 6.2566316036101626E-006 - 3.7163201120581581E-007 - 5.4155288561946149E-006 - 1.8786819174573623E-006 - 8.4967241033663055E-006 - 2.6961279691572182E-006 - 1.7416078313976980E-006 - 2.1226880367527847E-006 - 5.6985152337693010E-006 - 4.5338244070194161E-007 - 8.1224288329971305E-006 - 1.2583441000010966E-007 - 4.0920562284569457E-006 - 7.6329426691450179E-006 - 1.7160963129884001E-006 - 6.6076648110477140E-006 - 1.0731735278990562E-006 - 6.1073303601923391E-006 - 1.8860699966918250E-006 - 7.1975926015998006E-007 - 4.3303800633873823E-006 - 2.7581618400452237E-006 - 6.0028583379676688E-006 - 7.7027018498727923E-006 - 3.6999785309263525E-006 - 7.4314767294185678E-006 - 1.4088729129277634E-006 - 6.8721711355692943E-006 - 5.6081301435941245E-006 - 1.8073198064823138E-006 - 3.2030269920557565E-006 - 4.6370278037429165E-006 - 3.2322542793190805E-006 - 8.0124486124572907E-007 - 4.1995613708513501E-006 - 3.4756871649772262E-006 - 8.0509458007326399E-006 - 6.5522790234765574E-006 - 2.5617201839219564E-007 - 9.6199052717720368E-006 - 5.1882924569606104E-006 - 4.4263082100511957E-006 - 7.7019916315718217E-006 - 7.9953275316002850E-006 - 2.4508656586831879E-007 - 1.9660391811454640E-006 - 1.1556395663871900E-007 - 7.2925513202755867E-007 - 8.6334898352408556E-006 - 9.1764148916426310E-006 - 1.9877147326789627E-006 - 6.6911730452935485E-006 - 4.8461832878537443E-006 - 9.5998598209927511E-007 - 6.5297643919457324E-006 - 3.3382674828396801E-006 +-0.00008444 +-0.00003393 +0.00009230 +-0.00006224 +-0.00007552 +0.00009380 +0.00001181 +0.00009928 +-0.00007969 +0.00001493 +0.00007319 +-0.00004341 +-0.00002337 +0.00007838 +-0.00006033 +-0.00002512 +0.00004281 +0.00004631 +-0.00008111 +-0.00001758 +0.00002855 +0.00008134 +-0.00003731 +0.00002581 +0.00007723 +-0.00007482 +-0.00004745 +-0.00006915 +-0.00008698 +-0.00001448 +-0.00005086 +0.00001145 +-0.00003374 +0.00009208 +0.00007611 diff --git a/jastrow b/jastrow index 4279f94..51f8827 100755 Binary files a/jastrow and b/jastrow differ diff --git a/jastrow_provider.irp.f b/jastrow_provider.irp.f index 81a8635..9652d85 100644 --- a/jastrow_provider.irp.f +++ b/jastrow_provider.irp.f @@ -3,11 +3,12 @@ BEGIN_PROVIDER [ double precision, jastrow_full ] BEGIN_DOC ! Complete jastrow factor END_DOC + integer :: i, j print *, factor_ee print *, factor_en print *, factor_een - + jastrow_full = dexp(factor_ee + factor_en + factor_een) END_PROVIDER diff --git a/nuclei.irp.f b/nuclei.irp.f index 090132a..d11f19a 100644 --- a/nuclei.irp.f +++ b/nuclei.irp.f @@ -24,11 +24,13 @@ BEGIN_PROVIDER [ double precision, nuc_coord, (nnuc, 3) ] ! Nuclei coordinates END_DOC character(len=*), parameter :: FILE_NAME = "geometry.txt" - integer :: fu, rc + integer :: fu, rc, i open(action='read', file=FILE_NAME, iostat=rc, newunit=fu) - read(fu, *) nuc_coord + do i = 1, nnuc + read(fu, *) nuc_coord(i, :) + end do close(fu) @@ -42,12 +44,12 @@ BEGIN_PROVIDER [ double precision, elnuc_dist, (nelec, nnuc) ] integer :: i, j double precision :: x, y, z do j = 1, nnuc - do i = 1, nelec - x = elec_coord(i, 1) - nuc_coord(j, 1) - y = elec_coord(i, 2) - nuc_coord(j, 2) - z = elec_coord(i, 3) - nuc_coord(j, 3) - elnuc_dist(i, j) = dsqrt( x*x + y*y + z*z ) - enddo + do i = 1, nelec + x = elec_coord(i, 1) - nuc_coord(j, 1) + y = elec_coord(i, 2) - nuc_coord(j, 2) + z = elec_coord(i, 3) - nuc_coord(j, 3) + elnuc_dist(i, j) = dsqrt( x*x + y*y + z*z ) + enddo enddo END_PROVIDER @@ -64,10 +66,10 @@ BEGIN_PROVIDER [double precision, factor_en] do j = 1 , nnuc do i = 1, nelec - x = rescale_en(i, j) + x = rescale_en(i, j) * rescale_en(i, j) do p = 2, naord x = x * rescale_en(i, j) - pow_ser = pow_ser + aord_vect(p, typenuc_arr(j)) * x + pow_ser = pow_ser + aord_vect(p + 1, typenuc_arr(j)) * x end do factor_en = factor_en + aord_vect(1, typenuc_arr(j)) * rescale_en(i, j) & / (1 + aord_vect(2, typenuc_arr(j)) * rescale_en(i, j)) + pow_ser diff --git a/orders.irp.f b/orders.irp.f index de52116..df09f13 100644 --- a/orders.irp.f +++ b/orders.irp.f @@ -21,10 +21,35 @@ BEGIN_PROVIDER [integer, ncord] END_DOC ncord = 5 END_PROVIDER + +BEGIN_PROVIDER [integer, dim_cord_vect] + implicit none + BEGIN_DOC + ! Recomputes the length of the unique C coefficients + END_DOC + integer :: k, p, l, lmax -BEGIN_PROVIDER [double precision, aord_vect, (naord, typenuc)] -&BEGIN_PROVIDER [double precision, bord_vect, (nbord)] -&BEGIN_PROVIDER [double precision, cord_vect, (0:ncord , 0:ncord , ncord , typenuc)] + dim_cord_vect = 0 + + do p = 2, ncord + do k = 0, p - 1 + if ( k /= 0 ) then + lmax = p - k + else + lmax = p - k - 2 + end if + do l = iand(p - k, 1), lmax, 2 + dim_cord_vect = dim_cord_vect + 1 + end do + end do + end do + +END_PROVIDER + + +BEGIN_PROVIDER [double precision, aord_vect, (0:naord, typenuc)] +&BEGIN_PROVIDER [double precision, bord_vect, (0:nbord)] +&BEGIN_PROVIDER [double precision, cord_vect, (dim_cord_vect, typenuc)] implicit none BEGIN_DOC ! Read Jastow coefficients from file diff --git a/rescale.irp.f b/rescale.irp.f index ee91ec0..56b0d66 100644 --- a/rescale.irp.f +++ b/rescale.irp.f @@ -11,7 +11,7 @@ BEGIN_PROVIDER [ double precision, kappa_inv ] BEGIN_DOC ! inverse of kappa END_DOC - kappa_inv = 1.d0 / kappa + kappa_inv = 1.0d0 / kappa END_PROVIDER BEGIN_PROVIDER [ double precision, rescale_ee, (nelec, nelec) ] @@ -20,9 +20,10 @@ BEGIN_PROVIDER [ double precision, rescale_ee, (nelec, nelec) ] ! R = (1 - exp(-kappa r))/kappa for electron-electron for $J_{ee}$ END_DOC integer :: i, j - do j=1,nelec - do i=1,nelec - rescale_ee(i, j) = (1.d0 - dexp(-kappa * elec_dist(i, j))) * kappa_inv + + do j = 1, nelec + do i = 1, nelec + rescale_ee(i, j) = (1.0d0 - dexp(-kappa * elec_dist(i, j))) * kappa_inv enddo enddo END_PROVIDER @@ -33,6 +34,7 @@ BEGIN_PROVIDER [ double precision, rescale_en, (nelec, nnuc) ] ! R = (1 - exp(-kappa r))/kappa for electron-nucleus for $J_{en}$ END_DOC integer :: i, j + do j = 1, nnuc do i = 1, nelec rescale_en(i, j) = (1.d0 - dexp(-kappa * elnuc_dist(i, j))) * kappa_inv @@ -46,6 +48,7 @@ BEGIN_PROVIDER [double precision, rescale_een_e, (nelec, nelec)] ! R = exp(-kappa r) for electron-electron for $J_{een}$ END_DOC integer :: i, j + do j = 1, nelec do i = 1, nelec rescale_een_e(i, j) = dexp(-kappa * elec_dist(i, j)) @@ -59,6 +62,7 @@ BEGIN_PROVIDER [double precision, rescale_een_n, (nelec, nnuc)] ! R = exp(-kappa r) for electron-electron for $J_{een}$ END_DOC integer :: i, j + do j = 1, nnuc do i = 1, nelec rescale_een_n(i, j) = dexp(-kappa * elnuc_dist(i, j))