TREX/irpjast
1
0
Fork 0
mirror of https://github.com/TREX-CoE/irpjast.git synced 2024-07-23 03:07:47 +02:00
Code Issues Projects Releases Wiki Activity

Corrected distances

Browse Source
This commit is contained in:
Panadestein 2020-12-08 13:08:23 +01:00
parent 274b648e4c
commit ea0d11a899
9 changed files with 117 additions and 235 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
el_nuc_el.irp.f
View File

@ -3,7 +3,7 @@ BEGIN_PROVIDER [double precision, factor_een]
BEGIN_DOC BEGIN_DOC
! Electron-electron nucleus contribution to Jastrow factor ! Electron-electron nucleus contribution to Jastrow factor
END_DOC 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 double precision :: x, y, z, t, c_inv, u, a, b, a2, b2, c, t0
PROVIDE cord_vect PROVIDE cord_vect
@ -19,6 +19,7 @@ BEGIN_PROVIDER [double precision, factor_een]
b2 = b * b b2 = b * b
c = rescale_een_n(i, alpha) * rescale_een_n(j, alpha) c = rescale_een_n(i, alpha) * rescale_een_n(j, alpha)
c_inv = 1.0d0 / c c_inv = 1.0d0 / c
cindex = 0
do p = 2, ncord do p = 2, ncord
x = 1.0d0 x = 1.0d0
do k = 0, p - 1 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 ! This can be used in case of a flatten cord_vect
! cidx = 1 + l + (ncord + 1) * k + (ncord + 1) * (ncord + 1) * (p - 1) + & ! cidx = 1 + l + (ncord + 1) * k + (ncord + 1) * (ncord + 1) * (p - 1) + &
! (ncord + 1) * (ncord + 1) * ncord * (alpha - 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 t = t * c_inv
y = y * a2 y = y * a2
z = z * b2 z = z * b2

20
elec_coord.txt
View File

@ -1,10 +1,10 @@
0.47107089936082125 0.36934601008972345 0.13813129163530469 -0.250655104764153 0.503070975550133 -0.166554344502303
0.11734479025488309 0.36584274529092520 0.92672173271006275 -0.587812193472177 -0.128751981129274 0.187773606533075
0.35754797210216693 0.24703841579308339 0.17042939824115788 1.61335569047166 -0.615556732874863 -1.43165470979934
0.31813480208419698 0.81855922566341688 0.93785308019180225 -4.901239896295210E-003 -1.120440036458986E-002 1.99761909330422
0.69675380000569165 0.48537277058577866 2.8246759780116260E-002 0.766647499681200 -0.293515395797937 3.66454589201239
0.42646034508932040 5.1852406712776689E-002 0.93402091856761427 -0.127732483187947 -0.138975497694196 -8.669850480215846E-002
0.85418292218111624 0.51831660694819415 0.26738545106704503 -0.232271834949124 -1.059321673434182E-002 -0.504862241464867
0.34322893096420737 0.80132196716714288 0.33887619652274625 1.09360863531826 -2.036103063808752E-003 -2.702796910818986E-002
0.69754618540855695 0.39691148679374866 0.12716214781484358 -0.108090166832043 0.189161729653261 2.15398313919894
0.57090264858892203 0.70538842967387250 0.31300897404191264 0.397978144318712 -0.254277292595981 2.54553335476344

33
electrons.irp.f
View File

@ -13,14 +13,15 @@ BEGIN_PROVIDER [ double precision, elec_coord, (nelec, 3) ]
! Electron coordinates ! Electron coordinates
END_DOC END_DOC
character(len=*), parameter :: FILE_NAME = "elec_coord.txt" 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) 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) close(fu)
END_PROVIDER END_PROVIDER
BEGIN_PROVIDER [ double precision, elec_dist, (nelec, nelec) ] BEGIN_PROVIDER [ double precision, elec_dist, (nelec, nelec) ]
@ -28,15 +29,16 @@ BEGIN_PROVIDER [ double precision, elec_dist, (nelec, nelec) ]
BEGIN_DOC BEGIN_DOC
! e-e distance ! e-e distance
END_DOC END_DOC
integer :: i,j integer :: i, j
double precision :: x,y,z double precision :: x, y, z
do j = 1, nelec do j = 1, nelec
do i = 1, nelec do i = 1, nelec
x = elec_coord(i, 1) - elec_coord(j, 1) x = elec_coord(i, 1) - elec_coord(j, 1)
y = elec_coord(i, 2) - elec_coord(j, 2) y = elec_coord(i, 2) - elec_coord(j, 2)
z = elec_coord(i, 3) - elec_coord(j, 3) z = elec_coord(i, 3) - elec_coord(j, 3)
elec_dist(i, j) = dsqrt( x*x + y*y + z*z ) elec_dist(i, j) = dsqrt( x*x + y*y + z*z )
enddo enddo
enddo enddo
END_PROVIDER END_PROVIDER
@ -47,18 +49,19 @@ BEGIN_PROVIDER [double precision, factor_ee]
END_DOC END_DOC
integer :: i, j, p integer :: i, j, p
double precision :: pow_ser, x double precision :: pow_ser, x
factor_ee = 0.0d0 factor_ee = 0.0d0
pow_ser = 0.0d0 pow_ser = 0.0d0
do j = 1 , nelec do j = 1, nelec
do i = 1, nelec do i = 1, nelec
x = rescale_ee(i, j) x = rescale_ee(i, j) * rescale_ee(i, j)
do p = 2, nbord do p = 2, nbord
x = x * rescale_ee(i, j) 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 end do
factor_ee = factor_ee + bord_vect(1) * rescale_ee(i, j) & 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
end do end do

225
jast_coeffs.txt
View File

@ -1,190 +1,35 @@
4.7107089936082125E-004 -0.00008444
1.1734479025488309E-004 -0.00003393
3.5754797210216693E-004 0.00009230
3.1813480208419699E-004 -0.00006224
6.9675380000569162E-004 -0.00007552
4.7107089936082121E-008 0.00009380
1.1734479025488308E-008 0.00001181
3.5754797210216693E-008 0.00009928
3.1813480208419699E-008 -0.00007969
6.9675380000569162E-008 0.00001493
4.7107089936082132E-006 0.00007319
1.1734479025488310E-006 -0.00004341
3.5754797210216696E-006 -0.00002337
3.1813480208419699E-006 0.00007838
6.9675380000569168E-006 -0.00006033
4.2646034508932041E-006 -0.00002512
8.5418292218111638E-006 0.00004281
3.4322893096420741E-006 0.00004631
6.9754618540855701E-006 -0.00008111
5.7090264858892204E-006 -0.00001758
3.6934601008972349E-006 0.00002855
3.6584274529092522E-006 0.00008134
2.4703841579308343E-006 -0.00003731
8.1855922566341692E-006 0.00002581
4.8537277058577867E-006 0.00007723
5.1852406712776692E-007 -0.00007482
5.1831660694819417E-006 -0.00004745
8.0132196716714300E-006 -0.00006915
3.9691148679374872E-006 -0.00008698
7.0538842967387255E-006 -0.00001448
1.3813129163530470E-006 -0.00005086
9.2672173271006290E-006 0.00001145
1.7042939824115789E-006 -0.00003374
9.3785308019180231E-006 0.00009208
2.8246759780116260E-007 0.00007611
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

BIN
jastrow
View File

Binary file not shown.

1
jastrow_provider.irp.f
View File

@ -3,6 +3,7 @@ BEGIN_PROVIDER [ double precision, jastrow_full ]
BEGIN_DOC BEGIN_DOC
! Complete jastrow factor ! Complete jastrow factor
END_DOC END_DOC
integer :: i, j
print *, factor_ee print *, factor_ee
print *, factor_en print *, factor_en

22
nuclei.irp.f
View File

@ -24,11 +24,13 @@ BEGIN_PROVIDER [ double precision, nuc_coord, (nnuc, 3) ]
! Nuclei coordinates ! Nuclei coordinates
END_DOC END_DOC
character(len=*), parameter :: FILE_NAME = "geometry.txt" character(len=*), parameter :: FILE_NAME = "geometry.txt"
integer :: fu, rc integer :: fu, rc, i
open(action='read', file=FILE_NAME, iostat=rc, newunit=fu) 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) close(fu)
@ -42,12 +44,12 @@ BEGIN_PROVIDER [ double precision, elnuc_dist, (nelec, nnuc) ]
integer :: i, j integer :: i, j
double precision :: x, y, z double precision :: x, y, z
do j = 1, nnuc do j = 1, nnuc
do i = 1, nelec do i = 1, nelec
x = elec_coord(i, 1) - nuc_coord(j, 1) x = elec_coord(i, 1) - nuc_coord(j, 1)
y = elec_coord(i, 2) - nuc_coord(j, 2) y = elec_coord(i, 2) - nuc_coord(j, 2)
z = elec_coord(i, 3) - nuc_coord(j, 3) z = elec_coord(i, 3) - nuc_coord(j, 3)
elnuc_dist(i, j) = dsqrt( x*x + y*y + z*z ) elnuc_dist(i, j) = dsqrt( x*x + y*y + z*z )
enddo enddo
enddo enddo
END_PROVIDER END_PROVIDER
@ -64,10 +66,10 @@ BEGIN_PROVIDER [double precision, factor_en]
do j = 1 , nnuc do j = 1 , nnuc
do i = 1, nelec do i = 1, nelec
x = rescale_en(i, j) x = rescale_en(i, j) * rescale_en(i, j)
do p = 2, naord do p = 2, naord
x = x * rescale_en(i, j) 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 end do
factor_en = factor_en + aord_vect(1, typenuc_arr(j)) * rescale_en(i, j) & 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 / (1 + aord_vect(2, typenuc_arr(j)) * rescale_en(i, j)) + pow_ser

31
orders.irp.f
View File

@ -22,9 +22,34 @@ BEGIN_PROVIDER [integer, ncord]
ncord = 5 ncord = 5
END_PROVIDER END_PROVIDER
BEGIN_PROVIDER [double precision, aord_vect, (naord, typenuc)] BEGIN_PROVIDER [integer, dim_cord_vect]
&BEGIN_PROVIDER [double precision, bord_vect, (nbord)] implicit none
&BEGIN_PROVIDER [double precision, cord_vect, (0:ncord , 0:ncord , ncord , typenuc)] BEGIN_DOC
! Recomputes the length of the unique C coefficients
END_DOC
integer :: k, p, l, lmax
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 implicit none
BEGIN_DOC BEGIN_DOC
! Read Jastow coefficients from file ! Read Jastow coefficients from file

12
rescale.irp.f
View File

@ -11,7 +11,7 @@ BEGIN_PROVIDER [ double precision, kappa_inv ]
BEGIN_DOC BEGIN_DOC
! inverse of kappa ! inverse of kappa
END_DOC END_DOC
kappa_inv = 1.d0 / kappa kappa_inv = 1.0d0 / kappa
END_PROVIDER END_PROVIDER
BEGIN_PROVIDER [ double precision, rescale_ee, (nelec, nelec) ] 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}$ ! R = (1 - exp(-kappa r))/kappa for electron-electron for $J_{ee}$
END_DOC END_DOC
integer :: i, j integer :: i, j
do j=1,nelec
do i=1,nelec do j = 1, nelec
rescale_ee(i, j) = (1.d0 - dexp(-kappa * elec_dist(i, j))) * kappa_inv do i = 1, nelec
rescale_ee(i, j) = (1.0d0 - dexp(-kappa * elec_dist(i, j))) * kappa_inv
enddo enddo
enddo enddo
END_PROVIDER 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}$ ! R = (1 - exp(-kappa r))/kappa for electron-nucleus for $J_{en}$
END_DOC END_DOC
integer :: i, j integer :: i, j
do j = 1, nnuc do j = 1, nnuc
do i = 1, nelec do i = 1, nelec
rescale_en(i, j) = (1.d0 - dexp(-kappa * elnuc_dist(i, j))) * kappa_inv 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}$ ! R = exp(-kappa r) for electron-electron for $J_{een}$
END_DOC END_DOC
integer :: i, j integer :: i, j
do j = 1, nelec do j = 1, nelec
do i = 1, nelec do i = 1, nelec
rescale_een_e(i, j) = dexp(-kappa * elec_dist(i, j)) 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}$ ! R = exp(-kappa r) for electron-electron for $J_{een}$
END_DOC END_DOC
integer :: i, j integer :: i, j
do j = 1, nnuc do j = 1, nnuc
do i = 1, nelec do i = 1, nelec
rescale_een_n(i, j) = dexp(-kappa * elnuc_dist(i, j)) rescale_een_n(i, j) = dexp(-kappa * elnuc_dist(i, j))