TREX/irpjast
1
0
Fork 0
mirror of https://github.com/TREX-CoE/irpjast.git synced 2025-01-03 10:06:11 +01:00
Code Issues Projects Releases Wiki Activity

Merge pull request #1 from Panadestein/opt

Browse Source 
Opt
This commit is contained in:
Ramon L. PANADES-BARRUETA 2020-12-08 23:49:25 +01:00 committed by GitHub
commit 266e9dc655
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
12 changed files with 263 additions and 96 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
Makefile
View File

@ -1,8 +1,8 @@
IRPF90 = irpf90 -a -d
IRPF90 = irpf90 #-a -d
FC = gfortran
FCFLAGS= -O2 -ffree-line-length-none -I .
NINJA = ninja
AR = ar
ARCHIVE= ar crs
RANLIB = ranlib
SRC=

54
el_nuc_el.irp.f
View File

@ -3,26 +3,58 @@ 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 = 0
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
factor_een = 0.0d0
do alpha = 1, nnuc
do j = 1, nelec
b = rescale_een_n(j, alpha)
do i = 1, nelec
u = rescale_een_e(i, j)
a = rescale_een_n(i, alpha)
a2 = a * a
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
do k = p - 1, 0
if ( k == 0 ) then
lmax = p - k - 2
else
x = 1.0d0
do k = 0, p - 1
if ( k /= 0 ) then
lmax = p - k
else
lmax = p - k - 2
end if
do l = lmax, 0
if ( mod(p - k - l, 2) == 0 ) then
factor_een = factor_een + cord_vect(p, k, l) * rescale_een_e(i, j) ** k &
* (rescale_een_n(i, alpha) ** l + rescale_een_n(j, alpha) ** l) * &
(rescale_een_n(i, alpha) * rescale_een_n(j, alpha)) ** ((p - k - l) * 0.5d0)
end if
t = x
do l = 1, rshift(p - k, 1)
t = t * c
end do
! We have suppressed this from the following loop:
! if ( iand(p - k - l, 1) == 0 ) then
!
! Start from l=0 when p-k is even
! Start from l=1 when p-k is odd
if (iand(p - k, 1) == 0) then
y = 1.0d0
z = 1.0d0
else
y = a
z = b
endif
do l = iand(p - k, 1), lmax, 2
! 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)
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
end do
x = x * u
end do
end do
end do

10
elec_coord.txt Normal file
View File

@ -0,0 +1,10 @@
-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

50
electrons.irp.f
View File

@ -3,7 +3,15 @@ BEGIN_PROVIDER [ integer, nelec ]
BEGIN_DOC
! Number of electrons
END_DOC
nelec = 2
nelec = 10
END_PROVIDER
BEGIN_PROVIDER [ integer, nelec_up ]
implicit none
BEGIN_DOC
! Number of alpha and beta electrons
END_DOC
nelec_up = 5
END_PROVIDER
@ -12,13 +20,17 @@ BEGIN_PROVIDER [ double precision, elec_coord, (nelec, 3) ]
BEGIN_DOC
! Electron coordinates
END_DOC
integer :: i,j
do j = 1 , 3
character(len=*), parameter :: FILE_NAME = "elec_coord.txt"
integer :: fu, rc, i, j
open(action='read', file=FILE_NAME, iostat=rc, newunit=fu)
do i = 1, nelec
call random_number(elec_coord(i, j))
enddo
read(fu, *) elec_coord(i, :)
end do
close(fu)
END_PROVIDER
BEGIN_PROVIDER [ double precision, elec_dist, (nelec, nelec) ]
@ -28,6 +40,7 @@ BEGIN_PROVIDER [ double precision, elec_dist, (nelec, nelec) ]
END_DOC
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)
@ -43,18 +56,35 @@ BEGIN_PROVIDER [double precision, factor_ee]
BEGIN_DOC
! Electron-electron contribution to Jastrow factor
END_DOC
integer :: i, j
double precision :: pow_ser = 0.0d0
integer :: i, j, p, ipar
double precision :: pow_ser, x, spin_fact
factor_ee = 0.0d0
do j = 1, nelec
do i = 1, nelec
x = rescale_ee(i, j)
pow_ser = 0.0d0
spin_fact = 1.0d0
ipar = 0
do p = 2, nbord
pow_ser = pow_ser + bord_vect(p + 1) * rescale_ee(i, j) ** p
x = x * rescale_ee(i, j)
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
if ((i.le.nelec_up .and. j.le.nelec_up) .or. &
(i.gt.nelec_up .and. j.gt.nelec_up)) then
spin_fact = 0.5d0
ipar = 1
end if
factor_ee = factor_ee + spin_fact * bord_vect(1) * rescale_ee(i, j) &
/ (1.0d0 + bord_vect(2) * rescale_ee(i, j)) + pow_ser
end do
end do
factor_ee = 0.5d0 * factor_ee
END_PROVIDER

2
geometry.txt Normal file
View File

@ -0,0 +1,2 @@
0.000000 0.000000 0.000000
0.000000 0.000000 2.059801

35
jast_coeffs.txt Normal file
View File

@ -0,0 +1,35 @@
0.00000000
0.00000000
-0.380512
-0.157996
-0.031558
0.021512
0.5000000
0.153660
0.0672262
0.021570
0.0073096
0.002866
0.571702
-0.5142530
-0.513043
0.009486
-0.004205
0.4263258
0.0828815
0.0051186
-0.0029978
-0.0052704
-0.000075
-0.0830165
0.0145434
0.0514351
0.000925
-0.0040991
0.0043276
-0.00165447
0.002614
-0.001477
-0.0011370
-0.04010475
0.00610671

BIN
jastrow
View File

Binary file not shown.

2
jastrow.irp.f
View File

@ -1,6 +1,6 @@
program jastrow
implicit none
print *, 'The total Jastrow factor'
print *, dexp(factor_ee + factor_en + factor_een)
print *, jastrow_full
end program

15
jastrow_provider.irp.f Normal file
View File

@ -0,0 +1,15 @@
BEGIN_PROVIDER [ double precision, jastrow_full ]
implicit none
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

45
nuclei.irp.f
View File

@ -7,21 +7,36 @@ BEGIN_PROVIDER [ integer, nnuc ]
END_PROVIDER
BEGIN_PROVIDER [ integer, typenuc ]
&BEGIN_PROVIDER [integer, typenuc_arr, (nnuc)]
implicit none
BEGIN_DOC
! Number of nuclei
END_DOC
typenuc = 1
typenuc_arr = (/1, 1/)
END_PROVIDER
BEGIN_PROVIDER [ double precision, nuc_coord, (nnuc, 3) ]
implicit none
BEGIN_DOC
! Nuclei coordinates
END_DOC
integer :: i, j
do j = 1 , 3
character(len=*), parameter :: FILE_NAME = "geometry.txt"
integer :: fu, rc, i
open(action='read', file=FILE_NAME, iostat=rc, newunit=fu)
do i = 1, nnuc
call random_number(nuc_coord(i, j))
enddo
read(fu, *) nuc_coord(i, :)
end do
close(fu)
END_PROVIDER
BEGIN_PROVIDER [ double precision, elnuc_dist, (nnuc, nnuc) ]
BEGIN_PROVIDER [ double precision, elnuc_dist, (nelec, nnuc) ]
implicit none
BEGIN_DOC
! e-n distance
@ -43,19 +58,25 @@ BEGIN_PROVIDER [double precision, factor_en]
BEGIN_DOC
! Electron-nuclei contribution to Jastrow factor
END_DOC
integer :: i, j, p
double precision :: pow_ser = 0.0d0
integer :: i, j, p, q
double precision :: pow_ser, x
factor_en = 0.0d0
do j = 1 , nnuc
do i = 1, nnuc
do i = 1, nelec
x = rescale_en(i, j)
pow_ser = 0.0d0
do p = 2, naord
pow_ser = pow_ser + aord_vect(p + 1) * rescale_en(i, j) ** p
x = x * rescale_en(i, j)
pow_ser = pow_ser + aord_vect(p + 1, typenuc_arr(j)) * x
end do
factor_en = factor_en + aord_vect(1) * rescale_en(i, j) &
/ (1 + aord_vect(2) * rescale_en(i, j)) + pow_ser
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
end do
end do
factor_en = 0.5d0 * factor_en
END_PROVIDER

54
orders.irp.f
View File

@ -22,32 +22,50 @@ BEGIN_PROVIDER [integer, ncord]
ncord = 5
END_PROVIDER
BEGIN_PROVIDER [double precision, aord_vect, (naord)]
BEGIN_PROVIDER [integer, dim_cord_vect]
implicit none
BEGIN_DOC
! Vector of the `a' coefficients
! Recomputes the length of the unique C coefficients
END_DOC
do i = 1, naord
call random_number(aord_vect)
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, bord_vect, (nbord)]
implicit none
BEGIN_DOC
! Vector of the `b' coefficients
END_DOC
do i = 1, nbord
call random_number(bord_vect)
end do
END_PROVIDER
BEGIN_PROVIDER [double precision, cord_vect, (ncord)]
BEGIN_PROVIDER [double precision, aord_vect, (naord + 1, typenuc)]
&BEGIN_PROVIDER [double precision, bord_vect, (nbord + 1)]
&BEGIN_PROVIDER [double precision, cord_vect, (dim_cord_vect, typenuc)]
implicit none
BEGIN_DOC
! Vector of the `c' coefficients
! Read Jastow coefficients from file
END_DOC
do i = 1, ncord
call random_number(cord_vect)
end do
PROVIDE naord
PROVIDE nbord
PROVIDE ncord
character(len=*), parameter :: FILE_NAME = "jast_coeffs.txt"
integer :: i, fu, rc
open(action='read', file=FILE_NAME, iostat=rc, newunit=fu)
read(fu, *) aord_vect
read(fu, *) bord_vect
read(fu, *) cord_vect
close(fu)
END_PROVIDER

14
rescale.irp.f
View File

@ -3,7 +3,7 @@ BEGIN_PROVIDER [ double precision, kappa ]
BEGIN_DOC
! Constant in rescaling
END_DOC
kappa = 1.5d0
kappa = 0.6d0
END_PROVIDER
BEGIN_PROVIDER [ double precision, kappa_inv ]
@ -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
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
@ -40,12 +42,13 @@ BEGIN_PROVIDER [ double precision, rescale_en, (nelec, nnuc) ]
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, rescale_een_e, (nelec, 3)]
BEGIN_PROVIDER [double precision, rescale_een_e, (nelec, nelec)]
implicit none
BEGIN_DOC
! 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))
@ -53,12 +56,13 @@ BEGIN_PROVIDER [double precision, rescale_een_e, (nelec, 3)]
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, rescale_een_n, (nnuc, 3)]
BEGIN_PROVIDER [double precision, rescale_een_n, (nelec, nnuc)]
implicit none
BEGIN_DOC
! 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))