From bc39df1a0163ad49d721a980dbd8955cfc0c200a Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Tue, 18 Dec 2018 19:56:00 +0100 Subject: [PATCH] added density_for_dft, dft_keywords and data_energy_and_density --- src/ao_one_e_integrals/pot_ao_erf_ints.irp.f | 368 ++++++++++++++++++ src/becke_numerical_grid/README.rst | 9 +- .../{exemple.irp.f => example.irp.f} | 0 src/bitmask/bitmasks_routines.irp.f | 31 ++ src/bitmask/{exemple.irp.f => example.irp.f} | 0 src/data_energy_and_density/EZFIO.cfg | 26 ++ src/data_energy_and_density/NEED | 1 + src/data_energy_and_density/README.rst | 6 + .../save_one_body_dm.irp.f | 20 + src/density_for_dft/NEED | 3 + src/density_for_dft/README.rst | 4 + src/density_for_dft/density_for_dft.irp.f | 77 ++++ src/determinants/connected_to_ref.irp.f | 36 +- src/determinants/create_excitations.irp.f | 23 -- .../{exemple.irp.f => example.irp.f} | 0 src/determinants/excitations_utils.irp.f | 16 - src/dft_keywords/EZFIO.cfg | 37 ++ src/dft_keywords/NEED | 1 + src/dft_keywords/README.rst | 7 + src/dft_keywords/keywords.irp.f | 19 + 20 files changed, 638 insertions(+), 46 deletions(-) create mode 100644 src/ao_one_e_integrals/pot_ao_erf_ints.irp.f rename src/becke_numerical_grid/{exemple.irp.f => example.irp.f} (100%) rename src/bitmask/{exemple.irp.f => example.irp.f} (100%) create mode 100644 src/data_energy_and_density/EZFIO.cfg create mode 100644 src/data_energy_and_density/NEED create mode 100644 src/data_energy_and_density/README.rst create mode 100644 src/data_energy_and_density/save_one_body_dm.irp.f create mode 100644 src/density_for_dft/NEED create mode 100644 src/density_for_dft/README.rst create mode 100644 src/density_for_dft/density_for_dft.irp.f rename src/determinants/{exemple.irp.f => example.irp.f} (100%) delete mode 100644 src/determinants/excitations_utils.irp.f create mode 100644 src/dft_keywords/EZFIO.cfg create mode 100644 src/dft_keywords/NEED create mode 100644 src/dft_keywords/README.rst create mode 100644 src/dft_keywords/keywords.irp.f diff --git a/src/ao_one_e_integrals/pot_ao_erf_ints.irp.f b/src/ao_one_e_integrals/pot_ao_erf_ints.irp.f new file mode 100644 index 00000000..d076609c --- /dev/null +++ b/src/ao_one_e_integrals/pot_ao_erf_ints.irp.f @@ -0,0 +1,368 @@ +subroutine give_all_erf_kl_ao(integrals_ao,mu_in,C_center) + implicit none + BEGIN_DOC +! subroutine that returs all integrals over r of type erf(mu_in * |r-C_center|)/|r-C_center| + END_DOC + double precision, intent(in) :: mu_in,C_center(3) + double precision, intent(out) :: integrals_ao(ao_num,ao_num) + double precision :: NAI_pol_mult_erf_ao + integer :: i,j,l,k,m + do k = 1, ao_num + do m = 1, ao_num + integrals_ao(m,k) = NAI_pol_mult_erf_ao(m,k,mu_in,C_center) + enddo + enddo +end + + +double precision function NAI_pol_mult_erf_ao(i_ao,j_ao,mu_in,C_center) + implicit none + BEGIN_DOC + ! computes the following integral : + ! int[-infty;+infty] dr AO_i_ao (r) AO_j_ao(r) erf(mu_in * |r-C_center|)/|r-C_center| + END_DOC + integer, intent(in) :: i_ao,j_ao + double precision, intent(in) :: mu_in, C_center(3) + integer :: i,j,num_A,num_B, power_A(3), power_B(3), n_pt_in + double precision :: A_center(3), B_center(3),integral, alpha,beta, NAI_pol_mult_erf + num_A = ao_nucl(i_ao) + power_A(1:3)= ao_power(i_ao,1:3) + A_center(1:3) = nucl_coord(num_A,1:3) + num_B = ao_nucl(j_ao) + power_B(1:3)= ao_power(j_ao,1:3) + B_center(1:3) = nucl_coord(num_B,1:3) + n_pt_in = n_pt_max_integrals + NAI_pol_mult_erf_ao = 0.d0 + do i = 1, ao_prim_num(i_ao) + alpha = ao_expo_ordered_transp(i,i_ao) + do j = 1, ao_prim_num(j_ao) + beta = ao_expo_ordered_transp(j,j_ao) + integral = NAI_pol_mult_erf(A_center,B_center,power_A,power_B,alpha,beta,C_center,n_pt_in,mu_in) + NAI_pol_mult_erf_ao += integral * ao_coef_normalized_ordered_transp(j,j_ao)*ao_coef_normalized_ordered_transp(i,i_ao) + enddo + enddo + + + +end + + + +double precision function NAI_pol_mult_erf(A_center,B_center,power_A,power_B,alpha,beta,C_center,n_pt_in,mu_in) +! function that computes the folowing integral : +! int{dr} of (x-A_x)^ax (x-B_X)^bx exp(-alpha (x-A_x)^2 - beta (x-B_x)^2 ) erf(mu_in*(r-R_c))/(r-R_c) + +implicit none +integer, intent(in) :: n_pt_in +double precision,intent(in) :: C_center(3),A_center(3),B_center(3),alpha,beta,mu_in +integer, intent(in) :: power_A(3),power_B(3) +integer :: i,j,k,l,n_pt +double precision :: P_center(3) + +double precision :: d(0:n_pt_in),pouet,coeff,dist,const,pouet_2,factor +double precision :: I_n_special_exact,integrate_bourrin,I_n_bibi +double precision :: V_e_n,const_factor,dist_integral,tmp +double precision :: accu,rint,p_inv,p,rho,p_inv_2 +integer :: n_pt_out,lmax +include 'Utils/constants.include.F' + p = alpha + beta + p_inv = 1.d0/p + p_inv_2 = 0.5d0 * p_inv + rho = alpha * beta * p_inv + + dist = 0.d0 + dist_integral = 0.d0 + do i = 1, 3 + P_center(i) = (alpha * A_center(i) + beta * B_center(i)) * p_inv + dist += (A_center(i) - B_center(i))*(A_center(i) - B_center(i)) + dist_integral += (P_center(i) - C_center(i))*(P_center(i) - C_center(i)) + enddo + const_factor = dist*rho + if(const_factor > 80.d0)then + NAI_pol_mult_erf = 0.d0 + return + endif + double precision :: p_new + p_new = mu_in/dsqrt(p+ mu_in * mu_in) + factor = dexp(-const_factor) + coeff = dtwo_pi * factor * p_inv * p_new + lmax = 20 + + ! print*, "b" + do i = 0, n_pt_in + d(i) = 0.d0 + enddo + n_pt = 2 * ( (power_A(1) + power_B(1)) +(power_A(2) + power_B(2)) +(power_A(3) + power_B(3)) ) + const = p * dist_integral * p_new * p_new + if (n_pt == 0) then + pouet = rint(0,const) + NAI_pol_mult_erf = coeff * pouet + return + endif + +! call give_polynom_mult_center_mono_elec_erf(A_center,B_center,alpha,beta,power_A,power_B,C_center,n_pt_in,d,n_pt_out,mu_in) + p_new = p_new * p_new + call give_polynom_mult_center_mono_elec_erf_opt(A_center,B_center,alpha,beta,power_A,power_B,C_center,n_pt_in,d,n_pt_out,mu_in,p,p_inv,p_inv_2,p_new,P_center) + + + if(n_pt_out<0)then + NAI_pol_mult_erf = 0.d0 + return + endif + accu = 0.d0 + +! sum of integrals of type : int {t,[0,1]} exp-(rho.(P-Q)^2 * t^2) * t^i + do i =0 ,n_pt_out,2 + accu += d(i) * rint(i/2,const) + enddo + NAI_pol_mult_erf = accu * coeff + +end + + +subroutine give_polynom_mult_center_mono_elec_erf_opt(A_center,B_center,alpha,beta,power_A,power_B,C_center,n_pt_in,d,n_pt_out,mu_in,p,p_inv,p_inv_2,p_new,P_center) +!!!! subroutine that returns the explicit polynom in term of the "t" variable of the following polynomw :: +!!!! I_x1(a_x, d_x,p,q) * I_x1(a_y, d_y,p,q) * I_x1(a_z, d_z,p,q) +!!!! it is for the nuclear electron atraction +implicit none +integer, intent(in) :: n_pt_in +integer,intent(out) :: n_pt_out +double precision, intent(in) :: A_center(3), B_center(3),C_center(3),p,p_inv,p_inv_2,p_new,P_center(3) +double precision, intent(in) :: alpha,beta,mu_in +integer, intent(in) :: power_A(3), power_B(3) +integer :: a_x,b_x,a_y,b_y,a_z,b_z +double precision :: d(0:n_pt_in) +double precision :: d1(0:n_pt_in) +double precision :: d2(0:n_pt_in) +double precision :: d3(0:n_pt_in) +double precision :: accu + accu = 0.d0 +!COMPTEUR irp_rdtsc1 = irp_rdtsc() + ASSERT (n_pt_in > 1) + + double precision :: R1x(0:2), B01(0:2), R1xp(0:2),R2x(0:2) + R1x(0) = (P_center(1) - A_center(1)) + R1x(1) = 0.d0 + R1x(2) = -(P_center(1) - C_center(1))* p_new + ! R1x = (P_x - A_x) - (P_x - C_x) ( t * mu/sqrt(p+mu^2) )^2 + R1xp(0) = (P_center(1) - B_center(1)) + R1xp(1) = 0.d0 + R1xp(2) =-(P_center(1) - C_center(1))* p_new + !R1xp = (P_x - B_x) - (P_x - C_x) ( t * mu/sqrt(p+mu^2) )^2 + R2x(0) = p_inv_2 + R2x(1) = 0.d0 + R2x(2) = -p_inv_2* p_new + !R2x = 0.5 / p - 0.5/p ( t * mu/sqrt(p+mu^2) )^2 + do i = 0,n_pt_in + d(i) = 0.d0 + enddo + do i = 0,n_pt_in + d1(i) = 0.d0 + enddo + do i = 0,n_pt_in + d2(i) = 0.d0 + enddo + do i = 0,n_pt_in + d3(i) = 0.d0 + enddo + integer :: n_pt1,n_pt2,n_pt3,dim,i + n_pt1 = n_pt_in + n_pt2 = n_pt_in + n_pt3 = n_pt_in + a_x = power_A(1) + b_x = power_B(1) + call I_x1_pol_mult_mono_elec(a_x,b_x,R1x,R1xp,R2x,d1,n_pt1,n_pt_in) +! print*,'passed the first I_x1' + if(n_pt1<0)then + n_pt_out = -1 + do i = 0,n_pt_in + d(i) = 0.d0 + enddo + return + endif + + R1x(0) = (P_center(2) - A_center(2)) + R1x(1) = 0.d0 + R1x(2) = -(P_center(2) - C_center(2))* p_new + ! R1x = (P_x - A_x) - (P_x - C_x) ( t * mu/sqrt(p+mu^2) )^2 + R1xp(0) = (P_center(2) - B_center(2)) + R1xp(1) = 0.d0 + R1xp(2) =-(P_center(2) - C_center(2))* p_new + !R1xp = (P_x - B_x) - (P_x - C_x) ( t * mu/sqrt(p+mu^2) )^2 + a_y = power_A(2) + b_y = power_B(2) + call I_x1_pol_mult_mono_elec(a_y,b_y,R1x,R1xp,R2x,d2,n_pt2,n_pt_in) +! print*,'passed the second I_x1' + if(n_pt2<0)then + n_pt_out = -1 + do i = 0,n_pt_in + d(i) = 0.d0 + enddo + return + endif + + + R1x(0) = (P_center(3) - A_center(3)) + R1x(1) = 0.d0 + R1x(2) = -(P_center(3) - C_center(3))* p_new + ! R1x = (P_x - A_x) - (P_x - C_x) ( t * mu/sqrt(p+mu^2) )^2 + R1xp(0) = (P_center(3) - B_center(3)) + R1xp(1) = 0.d0 + R1xp(2) =-(P_center(3) - C_center(3))* p_new + !R2x = 0.5 / p - 0.5/p ( t * mu/sqrt(p+mu^2) )^2 + a_z = power_A(3) + b_z = power_B(3) + +! print*,'a_z = ',a_z +! print*,'b_z = ',b_z + call I_x1_pol_mult_mono_elec(a_z,b_z,R1x,R1xp,R2x,d3,n_pt3,n_pt_in) +! print*,'passed the third I_x1' + if(n_pt3<0)then + n_pt_out = -1 + do i = 0,n_pt_in + d(i) = 0.d0 + enddo + return + endif + integer :: n_pt_tmp + n_pt_tmp = 0 + call multiply_poly(d1,n_pt1,d2,n_pt2,d,n_pt_tmp) + do i = 0,n_pt_tmp + d1(i) = 0.d0 + enddo + n_pt_out = 0 + call multiply_poly(d ,n_pt_tmp ,d3,n_pt3,d1,n_pt_out) + do i = 0, n_pt_out + d(i) = d1(i) + enddo + +end + + + + +subroutine give_polynom_mult_center_mono_elec_erf(A_center,B_center,alpha,beta,power_A,power_B,C_center,n_pt_in,d,n_pt_out,mu_in) +!!!! subroutine that returns the explicit polynom in term of the "t" variable of the following polynomw :: +!!!! I_x1(a_x, d_x,p,q) * I_x1(a_y, d_y,p,q) * I_x1(a_z, d_z,p,q) +!!!! it is for the nuclear electron atraction +implicit none +integer, intent(in) :: n_pt_in +integer,intent(out) :: n_pt_out +double precision, intent(in) :: A_center(3), B_center(3),C_center(3) +double precision, intent(in) :: alpha,beta,mu_in +integer, intent(in) :: power_A(3), power_B(3) +integer :: a_x,b_x,a_y,b_y,a_z,b_z +double precision :: d(0:n_pt_in) +double precision :: d1(0:n_pt_in) +double precision :: d2(0:n_pt_in) +double precision :: d3(0:n_pt_in) +double precision :: accu, pq_inv, p10_1, p10_2, p01_1, p01_2 +double precision :: p,P_center(3),rho,p_inv,p_inv_2 + accu = 0.d0 +!COMPTEUR irp_rdtsc1 = irp_rdtsc() + ASSERT (n_pt_in > 1) + p = alpha+beta + p_inv = 1.d0/p + p_inv_2 = 0.5d0/p + do i =1, 3 + P_center(i) = (alpha * A_center(i) + beta * B_center(i)) * p_inv + enddo + + double precision :: R1x(0:2), B01(0:2), R1xp(0:2),R2x(0:2) + R1x(0) = (P_center(1) - A_center(1)) + R1x(1) = 0.d0 + R1x(2) = -(P_center(1) - C_center(1))* mu_in**2 / (p+mu_in*mu_in) + ! R1x = (P_x - A_x) - (P_x - C_x) ( t * mu/sqrt(p+mu^2) )^2 + R1xp(0) = (P_center(1) - B_center(1)) + R1xp(1) = 0.d0 + R1xp(2) =-(P_center(1) - C_center(1))* mu_in**2 / (p+mu_in*mu_in) + !R1xp = (P_x - B_x) - (P_x - C_x) ( t * mu/sqrt(p+mu^2) )^2 + R2x(0) = p_inv_2 + R2x(1) = 0.d0 + R2x(2) = -p_inv_2* mu_in**2 / (p+mu_in*mu_in) + !R2x = 0.5 / p - 0.5/p ( t * mu/sqrt(p+mu^2) )^2 + do i = 0,n_pt_in + d(i) = 0.d0 + enddo + do i = 0,n_pt_in + d1(i) = 0.d0 + enddo + do i = 0,n_pt_in + d2(i) = 0.d0 + enddo + do i = 0,n_pt_in + d3(i) = 0.d0 + enddo + integer :: n_pt1,n_pt2,n_pt3,dim,i + n_pt1 = n_pt_in + n_pt2 = n_pt_in + n_pt3 = n_pt_in + a_x = power_A(1) + b_x = power_B(1) + call I_x1_pol_mult_mono_elec(a_x,b_x,R1x,R1xp,R2x,d1,n_pt1,n_pt_in) +! print*,'passed the first I_x1' + if(n_pt1<0)then + n_pt_out = -1 + do i = 0,n_pt_in + d(i) = 0.d0 + enddo + return + endif + + R1x(0) = (P_center(2) - A_center(2)) + R1x(1) = 0.d0 + R1x(2) = -(P_center(2) - C_center(2))* mu_in**2 / (p+mu_in*mu_in) + ! R1x = (P_x - A_x) - (P_x - C_x) ( t * mu/sqrt(p+mu^2) )^2 + R1xp(0) = (P_center(2) - B_center(2)) + R1xp(1) = 0.d0 + R1xp(2) =-(P_center(2) - C_center(2))* mu_in**2 / (p+mu_in*mu_in) + !R1xp = (P_x - B_x) - (P_x - C_x) ( t * mu/sqrt(p+mu^2) )^2 + a_y = power_A(2) + b_y = power_B(2) + call I_x1_pol_mult_mono_elec(a_y,b_y,R1x,R1xp,R2x,d2,n_pt2,n_pt_in) +! print*,'passed the second I_x1' + if(n_pt2<0)then + n_pt_out = -1 + do i = 0,n_pt_in + d(i) = 0.d0 + enddo + return + endif + + + R1x(0) = (P_center(3) - A_center(3)) + R1x(1) = 0.d0 + R1x(2) = -(P_center(3) - C_center(3))* mu_in**2 / (p+mu_in*mu_in) + ! R1x = (P_x - A_x) - (P_x - C_x) ( t * mu/sqrt(p+mu^2) )^2 + R1xp(0) = (P_center(3) - B_center(3)) + R1xp(1) = 0.d0 + R1xp(2) =-(P_center(3) - C_center(3))* mu_in**2 / (p+mu_in*mu_in) + !R2x = 0.5 / p - 0.5/p ( t * mu/sqrt(p+mu^2) )^2 + a_z = power_A(3) + b_z = power_B(3) + +! print*,'a_z = ',a_z +! print*,'b_z = ',b_z + call I_x1_pol_mult_mono_elec(a_z,b_z,R1x,R1xp,R2x,d3,n_pt3,n_pt_in) +! print*,'passed the third I_x1' + if(n_pt3<0)then + n_pt_out = -1 + do i = 0,n_pt_in + d(i) = 0.d0 + enddo + return + endif + integer :: n_pt_tmp + n_pt_tmp = 0 + call multiply_poly(d1,n_pt1,d2,n_pt2,d,n_pt_tmp) + do i = 0,n_pt_tmp + d1(i) = 0.d0 + enddo + n_pt_out = 0 + call multiply_poly(d ,n_pt_tmp ,d3,n_pt3,d1,n_pt_out) + do i = 0, n_pt_out + d(i) = d1(i) + enddo + +end + diff --git a/src/becke_numerical_grid/README.rst b/src/becke_numerical_grid/README.rst index 43434514..a860eef4 100644 --- a/src/becke_numerical_grid/README.rst +++ b/src/becke_numerical_grid/README.rst @@ -2,7 +2,14 @@ Becke Numerical Grid ==================== -Lebedev-Laikov grids, using the code distributed through CCL (http://www.ccl.net/). +This module contains all quantities needed to build the Becke's grid used in general for DFT integration. Note that it can be used for whatever integration in R^3 as long as the functions to be integrated are mostly concentrated near the atomic regions. + +This grid is built as the reunion of a spherical grid around each atom. Each spherical grid contains a certain number of radial and angular points. + +For a simple example of how to use the grid, see example.irp.f. + +The spherical integration uses Lebedev-Laikov grids, which was used from the code distributed through CCL (http://www.ccl.net/). +See next section for explanations and citation policies. .. code-block:: text diff --git a/src/becke_numerical_grid/exemple.irp.f b/src/becke_numerical_grid/example.irp.f similarity index 100% rename from src/becke_numerical_grid/exemple.irp.f rename to src/becke_numerical_grid/example.irp.f diff --git a/src/bitmask/bitmasks_routines.irp.f b/src/bitmask/bitmasks_routines.irp.f index 34fb244c..97d8ec15 100644 --- a/src/bitmask/bitmasks_routines.irp.f +++ b/src/bitmask/bitmasks_routines.irp.f @@ -1,3 +1,34 @@ +subroutine set_bit_to_integer(i_physical,key,Nint) + use bitmasks + BEGIN_DOC + ! set to 1 the bit number i_physical in the bitstring key + END_DOC + implicit none + integer, intent(in) :: i_physical,Nint + integer(bit_kind), intent(inout) :: key(Nint) + integer :: k,j,i + k = ishft(i_physical-1,-bit_kind_shift)+1 + j = i_physical-ishft(k-1,bit_kind_shift)-1 + key(k) = ibset(key(k),j) +end + + +subroutine clear_bit_to_integer(i_physical,key,Nint) + use bitmasks + BEGIN_DOC + ! set to 0 the bit number i_physical in the bitstring key + END_DOC + implicit none + integer, intent(in) :: i_physical,Nint + integer(bit_kind), intent(inout) :: key(Nint) + integer :: k,j,i + k = ishft(i_physical-1,-bit_kind_shift)+1 + j = i_physical-ishft(k-1,bit_kind_shift)-1 + key(k) = ibclr(key(k),j) +end + + + subroutine bitstring_to_list( string, list, n_elements, Nint) use bitmasks implicit none diff --git a/src/bitmask/exemple.irp.f b/src/bitmask/example.irp.f similarity index 100% rename from src/bitmask/exemple.irp.f rename to src/bitmask/example.irp.f diff --git a/src/data_energy_and_density/EZFIO.cfg b/src/data_energy_and_density/EZFIO.cfg new file mode 100644 index 00000000..7126f257 --- /dev/null +++ b/src/data_energy_and_density/EZFIO.cfg @@ -0,0 +1,26 @@ +[data_energy_var] +type: double precision +doc: Variational energy computed with the wave function +interface: ezfio, provider +size: (determinants.n_states) + +[data_energy_proj] +type: double precision +doc: Projected energy computed with the wave function +interface: ezfio, provider +size: (determinants.n_states) + +[data_one_body_alpha_dm_mo] +interface: ezfio, provider +doc: Alpha one body density matrix on the MO basis computed with the wave function +type: double precision +size: (mo_basis.mo_tot_num,mo_basis.mo_tot_num,determinants.n_states) + + +[data_one_body_beta_dm_mo] +interface: ezfio, provider +doc: Beta one body density matrix on the MO basis computed with the wave function +type: double precision +size: (mo_basis.mo_tot_num,mo_basis.mo_tot_num,determinants.n_states) + + diff --git a/src/data_energy_and_density/NEED b/src/data_energy_and_density/NEED new file mode 100644 index 00000000..d3d4d2c7 --- /dev/null +++ b/src/data_energy_and_density/NEED @@ -0,0 +1 @@ +determinants diff --git a/src/data_energy_and_density/README.rst b/src/data_energy_and_density/README.rst new file mode 100644 index 00000000..8654090e --- /dev/null +++ b/src/data_energy_and_density/README.rst @@ -0,0 +1,6 @@ +============ +Data energy and density +============ + + +This module contains some global variables (such as densities and energies) which are stored in the EZFIO folder in a different place than determinants. This is used in practice to store density matrices which can be obtained from any methods, as long as they are stored in the same MO basis which is used for the calculations. In |RS-DFT| calculations, this can be done to perform damping on the density in order to speed up convergence. diff --git a/src/data_energy_and_density/save_one_body_dm.irp.f b/src/data_energy_and_density/save_one_body_dm.irp.f new file mode 100644 index 00000000..7d65e901 --- /dev/null +++ b/src/data_energy_and_density/save_one_body_dm.irp.f @@ -0,0 +1,20 @@ +program save_one_body_dm + implicit none + BEGIN_DOC +! programs that computes the one body density on the mo basis for alpha and beta electrons from the wave function stored in the EZFIO folder, and then save it into the EZFIO folder data_energy_and_density. +! +! Then, the global variable data_one_body_alpha_dm_mo and data_one_body_beta_dm_mo will automatically read the density in a further calculation. +! +! This can be used to perform dampin on the density in RS-DFT calculation (see the density_for_dft module). + END_DOC + read_wf = .True. + touch read_wf + call routine + +end + +subroutine routine + + call ezfio_set_data_energy_and_density_data_one_body_alpha_dm_mo(one_body_dm_mo_alpha) + call ezfio_set_data_energy_and_density_data_one_body_beta_dm_mo(one_body_dm_mo_beta) +end diff --git a/src/density_for_dft/NEED b/src/density_for_dft/NEED new file mode 100644 index 00000000..80fba7c9 --- /dev/null +++ b/src/density_for_dft/NEED @@ -0,0 +1,3 @@ +determinants +dft_keywords +data_energy_and_density diff --git a/src/density_for_dft/README.rst b/src/density_for_dft/README.rst new file mode 100644 index 00000000..4019f77f --- /dev/null +++ b/src/density_for_dft/README.rst @@ -0,0 +1,4 @@ +========== +DM_for_dft +========== + diff --git a/src/density_for_dft/density_for_dft.irp.f b/src/density_for_dft/density_for_dft.irp.f new file mode 100644 index 00000000..b106ca0b --- /dev/null +++ b/src/density_for_dft/density_for_dft.irp.f @@ -0,0 +1,77 @@ +BEGIN_PROVIDER [double precision, one_body_dm_mo_alpha_for_dft, (mo_tot_num,mo_tot_num, N_states)] + implicit none + BEGIN_DOC +! density matrix for alpha electrons in the MO basis used for all DFT calculations based on the density + END_DOC + double precision :: delta_alpha(mo_tot_num,mo_tot_num,N_states) + if(density_for_dft .EQ. "damping_rs_dft")then + delta_alpha = one_body_dm_mo_alpha - data_one_body_alpha_dm_mo + one_body_dm_mo_alpha_for_dft = data_one_body_alpha_dm_mo + damping_for_rs_dft * delta_alpha + else if (density_for_dft .EQ. "input_density")then + one_body_dm_mo_alpha_for_dft = data_one_body_alpha_dm_mo + else if (density_for_dft .EQ. "WFT")then + provide mo_coef + one_body_dm_mo_alpha_for_dft = one_body_dm_mo_alpha + endif + +END_PROVIDER + +BEGIN_PROVIDER [double precision, one_body_dm_mo_beta_for_dft, (mo_tot_num,mo_tot_num, N_states)] + implicit none + BEGIN_DOC +! density matrix for beta electrons in the MO basis used for all DFT calculations based on the density + END_DOC + double precision :: delta_beta(mo_tot_num,mo_tot_num,N_states) + if(density_for_dft .EQ. "damping_rs_dft")then + delta_beta = one_body_dm_mo_beta - data_one_body_beta_dm_mo + one_body_dm_mo_beta_for_dft = data_one_body_beta_dm_mo + damping_for_rs_dft * delta_beta + else if (density_for_dft .EQ. "input_density")then + one_body_dm_mo_beta_for_dft = data_one_body_beta_dm_mo + else if (density_for_dft .EQ. "WFT")then + provide mo_coef + one_body_dm_mo_beta_for_dft = one_body_dm_mo_beta + endif +END_PROVIDER + +BEGIN_PROVIDER [double precision, one_body_dm_mo_for_dft, (mo_tot_num,mo_tot_num, N_states)] + implicit none + one_body_dm_mo_for_dft = one_body_dm_mo_beta_for_dft + one_body_dm_mo_alpha_for_dft +END_PROVIDER + +BEGIN_PROVIDER [double precision, one_body_dm_average_mo_for_dft, (mo_tot_num,mo_tot_num)] + implicit none + integer :: i + one_body_dm_average_mo_for_dft = 0.d0 + do i = 1, N_states + one_body_dm_average_mo_for_dft(:,:) += one_body_dm_mo_for_dft(:,:,i) * state_average_weight(i) + enddo +END_PROVIDER + + BEGIN_PROVIDER [ double precision, one_body_dm_alpha_ao_for_dft, (ao_num,ao_num,N_states) ] +&BEGIN_PROVIDER [ double precision, one_body_dm_beta_ao_for_dft, (ao_num,ao_num,N_states) ] + BEGIN_DOC +! one body density matrix on the AO basis based on one_body_dm_mo_alpha_for_dft + END_DOC + implicit none + integer :: i,j,k,l,istate + double precision :: mo_alpha,mo_beta + + one_body_dm_alpha_ao_for_dft = 0.d0 + one_body_dm_beta_ao_for_dft = 0.d0 + do k = 1, ao_num + do l = 1, ao_num + do i = 1, mo_tot_num + do j = 1, mo_tot_num + do istate = 1, N_states + mo_alpha = one_body_dm_mo_alpha_for_dft(j,i,istate) + mo_beta = one_body_dm_mo_beta_for_dft(j,i,istate) + one_body_dm_alpha_ao_for_dft(l,k,istate) += mo_coef(k,i) * mo_coef(l,j) * mo_alpha + one_body_dm_beta_ao_for_dft(l,k,istate) += mo_coef(k,i) * mo_coef(l,j) * mo_beta + enddo + enddo + enddo + enddo + enddo + +END_PROVIDER + diff --git a/src/determinants/connected_to_ref.irp.f b/src/determinants/connected_to_ref.irp.f index 1b19ba47..07a94f9d 100644 --- a/src/determinants/connected_to_ref.irp.f +++ b/src/determinants/connected_to_ref.irp.f @@ -249,9 +249,21 @@ integer function connected_to_ref(key,keys,Nint,N_past_in,Ndet) logical :: t double precision :: hij_elec - ! output : 0 : not connected - ! i : connected to determinant i of the past - ! -i : is the ith determinant of the reference wf keys + BEGIN_DOC + ! input : key : a given Slater determinant + ! + ! : keys: a list of Slater determinants + ! + ! : Ndet: the number of Slater determinants in keys + ! + ! : N_past_in the number of Slater determinants for the connectivity research + ! + ! output : 0 : key not connected to the N_past_in first Slater determinants in keys + ! + ! i : key is connected to determinant i of keys + ! + ! -i : key is the ith determinant of the reference wf keys + END_DOC ASSERT (Nint > 0) ASSERT (Nint == N_int) @@ -349,9 +361,21 @@ integer function connected_to_ref_by_mono(key,keys,Nint,N_past_in,Ndet) logical :: t double precision :: hij_elec - ! output : 0 : not connected - ! i : connected to determinant i of the past - ! -i : is the ith determinant of the refernce wf keys + BEGIN_DOC + ! input : key : a given Slater determinant + ! + ! : keys: a list of Slater determinants + ! + ! : Ndet: the number of Slater determinants in keys + ! + ! : N_past_in the number of Slater determinants for the connectivity research + ! + ! output : 0 : key not connected by a MONO EXCITATION to the N_past_in first Slater determinants in keys + ! + ! i : key is connected by a MONO EXCITATION to determinant i of keys + ! + ! -i : key is the ith determinant of the reference wf keys + END_DOC ASSERT (Nint > 0) ASSERT (Nint == N_int) diff --git a/src/determinants/create_excitations.irp.f b/src/determinants/create_excitations.irp.f index fbcf0984..8551d44a 100644 --- a/src/determinants/create_excitations.irp.f +++ b/src/determinants/create_excitations.irp.f @@ -65,26 +65,3 @@ logical function is_spin_flip_possible(key_in,i_flip,ispin) endif end -subroutine set_bit_to_integer(i_physical,key,Nint) - use bitmasks - implicit none - integer, intent(in) :: i_physical,Nint - integer(bit_kind), intent(inout) :: key(Nint) - integer :: k,j,i - k = shiftr(i_physical-1,bit_kind_shift)+1 - j = i_physical-shiftl(k-1,bit_kind_shift)-1 - key(k) = ibset(key(k),j) -end - - -subroutine clear_bit_to_integer(i_physical,key,Nint) - use bitmasks - implicit none - integer, intent(in) :: i_physical,Nint - integer(bit_kind), intent(inout) :: key(Nint) - integer :: k,j,i - k = shiftr(i_physical-1,bit_kind_shift)+1 - j = i_physical-shiftl(k-1,bit_kind_shift)-1 - key(k) = ibclr(key(k),j) -end - diff --git a/src/determinants/exemple.irp.f b/src/determinants/example.irp.f similarity index 100% rename from src/determinants/exemple.irp.f rename to src/determinants/example.irp.f diff --git a/src/determinants/excitations_utils.irp.f b/src/determinants/excitations_utils.irp.f deleted file mode 100644 index a40b5690..00000000 --- a/src/determinants/excitations_utils.irp.f +++ /dev/null @@ -1,16 +0,0 @@ -subroutine apply_mono(i_hole,i_particle,ispin_excit,key_in,Nint) - implicit none - integer, intent(in) :: i_hole,i_particle,ispin_excit,Nint - integer(bit_kind), intent(inout) :: key_in(Nint,2) - integer :: k,j - use bitmasks - ! hole - k = shiftr(i_hole-1,bit_kind_shift)+1 - j = i_hole-shiftl(k-1,bit_kind_shift)-1 - key_in(k,ispin_excit) = ibclr(key_in(k,ispin_excit),j) - - k = shiftr(i_particle-1,bit_kind_shift)+1 - j = i_particle-shiftl(k-1,bit_kind_shift)-1 - key_in(k,ispin_excit) = ibset(key_in(k,ispin_excit),j) - -end diff --git a/src/dft_keywords/EZFIO.cfg b/src/dft_keywords/EZFIO.cfg new file mode 100644 index 00000000..c2eb92e5 --- /dev/null +++ b/src/dft_keywords/EZFIO.cfg @@ -0,0 +1,37 @@ +[exchange_functional] +type: character*(32) +doc: name of the exchange functional +interface: ezfio, provider, ocaml +default: short_range_LDA + +[correlation_functional] +type: character*(32) +doc: name of the correlation functional +interface: ezfio, provider, ocaml +default: short_range_LDA + +[HF_exchange] +type: double precision +doc: Percentage of HF exchange in the DFT model +interface: ezfio,provider,ocaml +default: 0. + +[mu_erf] +type: double precision +doc: cutting of the interaction in the range separated model +interface: ezfio,provider,ocaml +default: 0.5 +ezfio_name: mu_erf + + +[density_for_dft] +type: character*(32) +doc: Type of density used for DFT calculation. If WFT it uses the density of the WFT stored in terms of determinants. If input_density it uses the one-body dm stored in data_.../ . If damping_rs_dft it uses the damping density between WFT and input_density +interface: ezfio, provider, ocaml +default: WFT + +[damping_for_rs_dft] +type: double precision +doc: damping factor for the density used in RSFT. +interface: ezfio,provider,ocaml +default: 0.5 diff --git a/src/dft_keywords/NEED b/src/dft_keywords/NEED new file mode 100644 index 00000000..5a3182ed --- /dev/null +++ b/src/dft_keywords/NEED @@ -0,0 +1 @@ +ezfio_files diff --git a/src/dft_keywords/README.rst b/src/dft_keywords/README.rst new file mode 100644 index 00000000..089420cc --- /dev/null +++ b/src/dft_keywords/README.rst @@ -0,0 +1,7 @@ +============ +DFT Keywords +============ + +This module contains all keywords which are related to a DFT calculation + + diff --git a/src/dft_keywords/keywords.irp.f b/src/dft_keywords/keywords.irp.f new file mode 100644 index 00000000..10769c28 --- /dev/null +++ b/src/dft_keywords/keywords.irp.f @@ -0,0 +1,19 @@ +BEGIN_PROVIDER [ character*(32), DFT_TYPE] + implicit none + BEGIN_DOC +! defines the type of DFT applied: LDA, GGA etc ... + END_DOC + logical :: is_lda + if(correlation_functional.eq."None")then + is_lda = (index(exchange_functional,"LDA") .ne. 0) + else if(exchange_functional.eq."None")then + is_lda = (index(correlation_functional,"LDA") .ne. 0) + else + is_lda = (index(correlation_functional,"LDA") .ne. 0) .and. (index(exchange_functional,"LDA") .ne. 0) + endif + if(is_lda)then + DFT_TYPE = "LDA" + else + DFT_TYPE = "GGA" + endif +END_PROVIDER