9
1
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-12-22 03:23:29 +01:00
This commit is contained in:
Anthony Scemama 2020-04-28 15:49:18 +02:00
commit ef9e61bd24
9 changed files with 90 additions and 10 deletions

View File

@ -11,9 +11,6 @@
ELSE
ao_one_e_integrals = ao_integrals_n_e + ao_kinetic_integrals
IF (DO_PSEUDO) THEN
ao_one_e_integrals += ao_pseudo_integrals
ENDIF
ENDIF
DO j = 1, ao_num

View File

@ -76,6 +76,12 @@ BEGIN_PROVIDER [ double precision, ao_integrals_n_e, (ao_num,ao_num)]
!$OMP END DO
!$OMP END PARALLEL
endif
IF (DO_PSEUDO) THEN
ao_integrals_n_e += ao_pseudo_integrals
ENDIF
if (write_ao_integrals_n_e) then
call ezfio_set_ao_one_e_ints_ao_integrals_n_e(ao_integrals_n_e)
print *, 'AO N-e integrals written to disk'

View File

@ -26,6 +26,14 @@ subroutine ec_md_pbe_on_top_general(mu,rho_a,rho_b,grad_rho_a,grad_rho_b,on_top,
pi = 4.d0 * datan(1.d0)
eps_c_md_on_top_PBE = 0.d0
! convertion from (alpha,beta) formalism to (closed, open) formalism for the density
call rho_ab_to_rho_oc(rho_a,rho_b,rhoo,rhoc)
if(rhoc.lt.1.d-10)then
return
else if(on_top/(rhoc**2) .lt. 1.d-6)then
return
endif
grad_rho_a_2 = 0.d0
grad_rho_b_2 = 0.d0
grad_rho_a_b = 0.d0
@ -34,8 +42,7 @@ subroutine ec_md_pbe_on_top_general(mu,rho_a,rho_b,grad_rho_a,grad_rho_b,on_top,
grad_rho_b_2 += grad_rho_b(m)*grad_rho_b(m)
grad_rho_a_b += grad_rho_a(m)*grad_rho_b(m)
enddo
! convertion from (alpha,beta) formalism to (closed, open) formalism
call rho_ab_to_rho_oc(rho_a,rho_b,rhoo,rhoc)
! same same for gradients : convertion from (alpha,beta) formalism to (closed, open) formalism
call grad_rho_ab_to_grad_rho_oc(grad_rho_a_2,grad_rho_b_2,grad_rho_a_b,sigmaoo,sigmacc,sigmaco)
! usual PBE correlation energy using the density, spin polarization and density gradients for alpha/beta electrons

View File

@ -12,10 +12,6 @@ BEGIN_PROVIDER [ double precision, mo_one_e_integrals,(mo_num,mo_num)]
ELSE
mo_one_e_integrals = mo_integrals_n_e + mo_kinetic_integrals
IF (DO_PSEUDO) THEN
mo_one_e_integrals += mo_pseudo_integrals
ENDIF
ENDIF
IF (write_mo_one_e_integrals) THEN

View File

@ -311,3 +311,34 @@ BEGIN_PROVIDER [double precision, two_e_int_ii_f, (n_basis_orb,n_basis_orb,n_ina
enddo
END_PROVIDER
subroutine give_mu_of_r_cas(r,istate,mu_of_r,f_psi,n2_psi)
implicit none
BEGIN_DOC
! returns mu(r), f_psi, n2_psi for a general cas wave function
END_DOC
integer, intent(in) :: istate
double precision, intent(in) :: r(3)
double precision, intent(out) :: mu_of_r,f_psi,n2_psi
double precision :: f_ii_val_ab,two_bod_dens_ii
double precision :: f_ia_val_ab,two_bod_dens_ia
double precision :: f_aa_val_ab,two_bod_dens_aa
double precision :: sqpi,w_psi
sqpi = dsqrt(dacos(-1.d0))
! inactive-inactive part of f_psi(r1,r2)
call give_f_ii_val_ab(r,r,f_ii_val_ab,two_bod_dens_ii)
! inactive-active part of f_psi(r1,r2)
call give_f_ia_val_ab(r,r,f_ia_val_ab,two_bod_dens_ia,istate)
! active-active part of f_psi(r1,r2)
call give_f_aa_val_ab(r,r,f_aa_val_ab,two_bod_dens_aa,istate)
f_psi = f_ii_val_ab + f_ia_val_ab + f_aa_val_ab
n2_psi = two_bod_dens_ii + two_bod_dens_ia + two_bod_dens_aa
if(n2_psi.le.1.d-12.or.f_psi.le.0.d0.or.f_psi * n2_psi.lt.0.d0)then
w_psi = 1.d+10
else
w_psi = f_psi / n2_psi
endif
mu_of_r = w_psi * sqpi * 0.5d0
end

View File

@ -25,6 +25,7 @@
print*,'Providing act_2_rdm_ab_mo '
ispin = 3
act_2_rdm_ab_mo = 0.d0
provide mo_two_e_integrals_in_map
call wall_time(wall_1)
if(read_two_body_rdm_ab)then
print*,'Reading act_2_rdm_ab_mo from disk ...'

View File

@ -19,6 +19,7 @@ subroutine orb_range_2_rdm_openmp(big_array,dim1,norb,list_orb,ispin,u_0,N_st,sz
integer :: k
double precision, allocatable :: u_t(:,:)
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: u_t
PROVIDE mo_two_e_integrals_in_map
allocate(u_t(N_st,N_det))
do k=1,N_st
call dset_order(u_0(1,k),psi_bilinear_matrix_order,N_det)

41
src/utils/prim_in_r.irp.f Normal file
View File

@ -0,0 +1,41 @@
double precision function primitive_value_explicit(power_prim,center_prim,alpha,r)
implicit none
BEGIN_DOC
! Evaluates at "r" a primitive of type :
! (x - center_prim(1))**power_prim(1) (y - center_prim(2))**power_prim(2) * (z - center_prim(3))**power_prim(3)
!
! exp(-alpha * [(x - center_prim(1))**2 + (y - center_prim(2))**2 + (z - center_prim(3))**2] )
END_DOC
integer, intent(in) :: power_prim(3)
double precision, intent(in) :: center_prim(3),alpha
double precision, intent(in) :: r(3)
double precision :: dx,dy,dz,r2
dx = (r(1) - center_prim(1))
dy = (r(2) - center_prim(2))
dz = (r(3) - center_prim(3))
r2 = dx*dx + dy*dy + dz*dz
dx = dx**power_prim(1)
dy = dy**power_prim(2)
dz = dz**power_prim(3)
primitive_value_explicit = dexp(-alpha*r2) * dx * dy * dz
end
double precision function give_pol_in_r(r,pol,center, alpha,iorder, max_dim)
double precision :: r(3), center(3), alpha,pol(0:max_dim,3)
integer, intent(in) :: iorder(3), max_dim
integer :: i,m
double precision :: gauss(3), x
gauss = 0.d0
do m = 1, 3
x = r(m) - center(m)
do i = 0, iorder(m)
gauss(m) += pol(i,m) * dexp(-alpha *x**2 ) * x**i
enddo
enddo
give_pol_in_r = gauss(1) * gauss(2) * gauss(3)
end

View File

@ -589,11 +589,11 @@ subroutine end_parallel_job(zmq_to_qp_run_socket,zmq_socket_pull,name_in)
rc = f77_zmq_send(zmq_to_qp_run_socket, 'end_job '//trim(zmq_state),8+len(trim(zmq_state)),0)
rc = f77_zmq_recv(zmq_to_qp_run_socket, message, 512, 0)
if (trim(message(1:13)) == 'error waiting') then
call sleep(1)
cycle
else if (message(1:2) == 'ok') then
exit
endif
call sleep(1)
end do
if (i==0) then
print *, '.. Forcing kill ..'