jast 201 added

This commit is contained in:
AbdAmmar 2023-05-06 18:25:06 +02:00
parent 5eb3b69873
commit 868d5a162b
5 changed files with 120 additions and 10 deletions

View File

@ -9,6 +9,12 @@ doc: Coefficient of the i-th |AO| on the j-th |MO|
interface: ezfio
size: (ao_basis.ao_num,mo_basis.mo_num)
[mo_coef_aux]
type: double precision
doc: AUX Coefficient of the i-th |AO| on the j-th |MO|
interface: ezfio
size: (ao_basis.ao_num,mo_basis.mo_num)
[mo_coef_imag]
type: double precision
doc: Imaginary part of the MO coefficient of the i-th |AO| on the j-th |MO|

View File

@ -0,0 +1,53 @@
! ---
BEGIN_PROVIDER [double precision, mo_coef_aux, (ao_num,mo_num)]
implicit none
integer :: i, j
logical :: exists
double precision, allocatable :: buffer(:,:)
PROVIDE ezfio_filename
if (mpi_master) then
! Coefs
call ezfio_has_mo_basis_mo_coef_aux(exists)
endif
IRP_IF MPI_DEBUG
print *, irp_here, mpi_rank
call MPI_BARRIER(MPI_COMM_WORLD, ierr)
IRP_ENDIF
IRP_IF MPI
include 'mpif.h'
integer :: ierr
call MPI_BCAST(exists, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr)
if (ierr /= MPI_SUCCESS) then
stop 'Unable to read mo_coef_aux with MPI'
endif
IRP_ENDIF
if (exists) then
if (mpi_master) then
call ezfio_get_mo_basis_mo_coef_aux(mo_coef_aux)
write(*,*) 'Read mo_coef_aux'
endif
IRP_IF MPI
call MPI_BCAST(mo_coef_aux, mo_num*ao_num, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
if (ierr /= MPI_SUCCESS) then
stop 'Unable to read mo_coef_aux with MPI'
endif
IRP_ENDIF
else
! Orthonormalized AO basis
do i = 1, mo_num
do j = 1, ao_num
mo_coef_aux(j,i) = ao_ortho_canonical_coef(j,i)
enddo
enddo
endif
END_PROVIDER

View File

@ -408,24 +408,71 @@ subroutine mu_r_val_and_grad(r1, r2, mu_val, mu_der)
implicit none
double precision, intent(in) :: r1(3), r2(3)
double precision, intent(out) :: mu_val, mu_der(3)
double precision :: r(3)
double precision :: dm_a(1), dm_b(1), grad_dm_a(3,1), grad_dm_b(3,1)
double precision :: dm_tot, tmp1, tmp2, tmp3
if(j1b_type .eq. 200) then
double precision :: r(3), dm_a(1), dm_b(1), grad_dm_a(3,1), grad_dm_b(3,1)
double precision :: dm_tot, tmp
!
! r = 0.5 (r1 + r2)
!
! mu[rho(r)] = alpha sqrt(rho) + mu0 exp(-rho)
!
! d mu[rho(r)] / dx1 = 0.5 d mu[rho(r)] / dx
! d mu[rho(r)] / dx = [0.5 alpha / sqrt(rho) - mu0 exp(-rho)] (d rho(r) / dx)
!
PROVIDE mu_r_ct mu_erf
PROVIDE mu_r_ct
r(1) = 0.5d0 * (r1(1) + r2(1))
r(2) = 0.5d0 * (r1(2) + r2(2))
r(3) = 0.5d0 * (r1(3) + r2(3))
call density_and_grad_alpha_beta(r, dm_a, dm_b, grad_dm_a, grad_dm_b)
dm_tot = dm_a(1) + dm_b(1)
mu_val = mu_r_ct * dsqrt(dm_tot)
tmp = 0.25d0 * mu_r_ct / dm_tot
mu_der(1) = tmp * (grad_dm_a(1,1) + grad_dm_b(1,1))
mu_der(2) = tmp * (grad_dm_a(2,1) + grad_dm_b(2,1))
mu_der(3) = tmp * (grad_dm_a(3,1) + grad_dm_b(3,1))
dm_tot = dm_a(1) + dm_b(1)
tmp1 = dsqrt(dm_tot)
tmp2 = mu_erf * dexp(-dm_tot)
mu_val = mu_r_ct * tmp1 + tmp2
mu_der = 0.d0
if(dm_tot .lt. 1d-7) return
tmp3 = 0.25d0 * mu_r_ct / tmp1 - 0.5d0 * tmp2
mu_der(1) = tmp3 * (grad_dm_a(1,1) + grad_dm_b(1,1))
mu_der(2) = tmp3 * (grad_dm_a(2,1) + grad_dm_b(2,1))
mu_der(3) = tmp3 * (grad_dm_a(3,1) + grad_dm_b(3,1))
elseif(j1b_type .eq. 201) then
!
! r = 0.5 (r1 + r2)
!
! mu[rho(r)] = alpha rho + mu0 exp(-rho)
!
! d mu[rho(r)] / dx1 = 0.5 d mu[rho(r)] / dx
! d mu[rho(r)] / dx = [0.5 alpha / sqrt(rho) - mu0 exp(-rho)] (d rho(r) / dx)
!
PROVIDE mu_r_ct mu_erf
r(1) = 0.5d0 * (r1(1) + r2(1))
r(2) = 0.5d0 * (r1(2) + r2(2))
r(3) = 0.5d0 * (r1(3) + r2(3))
call density_and_grad_alpha_beta(r, dm_a, dm_b, grad_dm_a, grad_dm_b)
dm_tot = dm_a(1) + dm_b(1)
tmp2 = mu_erf * dexp(-dm_tot)
mu_val = mu_r_ct * dm_tot + tmp2
tmp3 = 0.5d0 * (mu_r_ct - tmp2)
mu_der(1) = tmp3 * (grad_dm_a(1,1) + grad_dm_b(1,1))
mu_der(2) = tmp3 * (grad_dm_a(2,1) + grad_dm_b(2,1))
mu_der(3) = tmp3 * (grad_dm_a(3,1) + grad_dm_b(3,1))
else

View File

@ -9,6 +9,10 @@ program print_tc_energy
my_n_pt_a_grid = 50
read_wf = .True.
touch read_wf
PROVIDE j1b_type
print*, 'j1b_type = ', j1b_type
touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
call write_tc_energy
end

View File

@ -176,7 +176,7 @@ default: 1.e-7
type: logical
doc: If |true|, the integrals of the three-body jastrow are computed with cycles
interface: ezfio,provider,ocaml
default: True
default: Flase
[thresh_biorthog_diag]
type: Threshold