9
1
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-12-22 19:43:32 +01:00

Removed dupicate provider

This commit is contained in:
Anthony Scemama 2024-12-17 10:47:31 +01:00
parent b3f481b11a
commit a61f6291da
26 changed files with 1412 additions and 10 deletions

View File

@ -0,0 +1,2 @@
Transcorrelated SCF
===================

View File

@ -0,0 +1,104 @@
[ao_extra_basis]
type: character*(256)
doc: Name of the |ao_extra| basis set
interface: ezfio
[ao_extra_only_1s]
type: logical
doc: If |true|, you know that the additional AO basis is built only with 1s functions
interface: ezfio, provider
default: true
[ao_extra_num]
type: integer
doc: Number of |ao_extras|
interface: ezfio, provider
[ao_extra_prim_num]
type: integer
doc: Number of primitives per |ao_extra|
size: (ao_extra_basis.ao_extra_num)
interface: ezfio, provider
[ao_extra_prim_num_max]
type: integer
doc: Maximum number of primitives
default: =maxval(ao_extra_basis.ao_extra_prim_num)
interface: ezfio
[ao_extra_nucl]
type: integer
doc: Index of the nucleus on which the |ao_extra| is centered
size: (ao_extra_basis.ao_extra_num)
interface: ezfio, provider
[ao_extra_power]
type: integer
doc: Powers of x, y and z for each |ao_extra|
size: (ao_extra_basis.ao_extra_num,3)
interface: ezfio, provider
[ao_extra_coef]
type: double precision
doc: Primitive coefficients, read from input. Those should not be used directly, as the MOs are expressed on the basis of **normalized** ao_extras.
size: (ao_extra_basis.ao_extra_num,ao_extra_basis.ao_extra_prim_num_max)
interface: ezfio, provider
[ao_extra_expo]
type: double precision
doc: Exponents for each primitive of each |ao_extra|
size: (ao_extra_basis.ao_extra_num,ao_extra_basis.ao_extra_prim_num_max)
interface: ezfio, provider
[ao_extra_md5]
type: character*(32)
doc: MD5 key, specific of the |ao_extra| basis
interface: ezfio, provider
[ao_extra_cartesian]
type: logical
doc: If |true|, use |ao_extras| in Cartesian coordinates (6d,10f,...)
interface: ezfio, provider
default: false
[ao_extra_normalized]
type: logical
doc: Use normalized basis functions
interface: ezfio, provider
default: true
[primitives_normalized_extra]
type: logical
doc: Use normalized primitive functions
interface: ezfio, provider
default: true
[ao_extra_expo_im]
type: double precision
doc: imag part for Exponents for each primitive of each cGTOs |ao_extra|
size: (ao_extra_basis.ao_extra_num,ao_extra_basis.ao_extra_prim_num_max)
interface: ezfio, provider
[ao_extra_expo_pw]
type: double precision
doc: plane wave part for each primitive GTOs |ao_extra|
size: (3,ao_extra_basis.ao_extra_num,ao_extra_basis.ao_extra_prim_num_max)
interface: ezfio, provider
[ao_extra_expo_phase]
type: double precision
doc: phase shift for each primitive GTOs |ao_extra|
size: (3,ao_extra_basis.ao_extra_num,ao_extra_basis.ao_extra_prim_num_max)
interface: ezfio, provider
[ao_extra_one_e_dm]
type: double precision
doc: reduced density matrix on the ao extra basis
size: (ao_extra_basis.ao_extra_num,ao_extra_basis.ao_extra_num)
interface: ezfio, provider
[ao_extra_center]
type: double precision
doc: shift with which the atoms are shifted to mimick p functions
interface: ezfio

View File

@ -0,0 +1,4 @@
2
H 0. 0. 0.
Li 0. 0. 1.0

3
src/ao_extra_basis/NEED Normal file
View File

@ -0,0 +1,3 @@
extra_nuclei
basis
ao_basis

View File

@ -0,0 +1,15 @@
===========
extra_basis
===========
Plugin to handle an extra basis, which is attached to the extra_nuclei.
It is essentially a duplication of all important quantities (coefficients, exponents and so on) of the usual |AO| basis.
An interesting feature is the possibility to fit any basis made at most with "p" functions onto a purely "s" basis.
This is done with the various scripts here:
- qp_fit_1s_basis : script that creates an |EZFIO| folder corresponding to an .xyz file and a basis fitted with only "s" functions
- qp_add_extra_fit_system : script that takes as input an |EZFIO| folder and an .xyz file and add an extra_basis and extra_nuclei with a 1s fitted basis
Ex:
qp_add_extra_fit_system LiH.ezfio/ h2o.xyz # takes the EZFIO folder "LiH.ezfio" and creates all necessary additional basis and nuclei based on h2o.xyz, but only with 1s functions.

View File

@ -0,0 +1,325 @@
BEGIN_PROVIDER [ integer, ao_extra_prim_num_max ]
implicit none
BEGIN_DOC
! Max number of primitives.
END_DOC
ao_extra_prim_num_max = maxval(ao_extra_prim_num)
END_PROVIDER
BEGIN_PROVIDER [ integer, ao_extra_shell, (ao_extra_num) ]
implicit none
BEGIN_DOC
! Index of the shell to which the ao_extra corresponds
END_DOC
integer :: i, j, k, n
k=0
do i=1,shell_num
n = shell_ang_mom(i)+1
do j=1,(n*(n+1))/2
k = k+1
ao_extra_shell(k) = i
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ integer, ao_extra_first_of_shell, (shell_num) ]
implicit none
BEGIN_DOC
! Index of the shell to which the ao_extra corresponds
END_DOC
integer :: i, j, k, n
k=1
do i=1,shell_num
ao_extra_first_of_shell(i) = k
n = shell_ang_mom(i)+1
k = k+(n*(n+1))/2
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_extra_coef_normalized, (ao_extra_num,ao_extra_prim_num_max) ]
&BEGIN_PROVIDER [ double precision, ao_extra_coef_normalization_factor, (ao_extra_num) ]
implicit none
BEGIN_DOC
! Coefficients including the |ao_extra| normalization
END_DOC
double precision :: norm,overlap_x,overlap_y,overlap_z,C_A(3), c
integer :: l, powA(3), nz
integer :: i,j,k
nz=100
C_A(1) = 0.d0
C_A(2) = 0.d0
C_A(3) = 0.d0
ao_extra_coef_normalized = 0.d0
do i=1,ao_extra_num
! powA(1) = ao_extra_power(i,1) + ao_extra_power(i,2) + ao_extra_power(i,3)
! powA(2) = 0
! powA(3) = 0
powA(1) = ao_extra_power(i,1)
powA(2) = ao_extra_power(i,2)
powA(3) = ao_extra_power(i,3)
! Normalization of the primitives
if (primitives_normalized_extra) then
do j=1,ao_extra_prim_num(i)
call overlap_gaussian_xyz(C_A,C_A,ao_extra_expo(i,j),ao_extra_expo(i,j), &
powA,powA,overlap_x,overlap_y,overlap_z,norm,nz)
ao_extra_coef_normalized(i,j) = ao_extra_coef(i,j)/dsqrt(norm)
enddo
else
do j=1,ao_extra_prim_num(i)
ao_extra_coef_normalized(i,j) = ao_extra_coef(i,j)
enddo
endif
! Normalization of the contracted basis functions
if (ao_extra_normalized) then
norm = 0.d0
do j=1,ao_extra_prim_num(i)
do k=1,ao_extra_prim_num(i)
call overlap_gaussian_xyz(C_A,C_A,ao_extra_expo(i,j),ao_extra_expo(i,k),powA,powA,overlap_x,overlap_y,overlap_z,c,nz)
norm = norm+c*ao_extra_coef_normalized(i,j)*ao_extra_coef_normalized(i,k)
enddo
enddo
ao_extra_coef_normalization_factor(i) = 1.d0/dsqrt(norm)
else
ao_extra_coef_normalization_factor(i) = 1.d0
endif
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_extra_coef_normalized_ordered, (ao_extra_num,ao_extra_prim_num_max) ]
&BEGIN_PROVIDER [ double precision, ao_extra_expo_ordered, (ao_extra_num,ao_extra_prim_num_max) ]
implicit none
BEGIN_DOC
! Sorted primitives to accelerate 4 index |MO| transformation
END_DOC
integer :: iorder(ao_extra_prim_num_max)
double precision :: d(ao_extra_prim_num_max,2)
integer :: i,j
do i=1,ao_extra_num
do j=1,ao_extra_prim_num(i)
iorder(j) = j
d(j,1) = ao_extra_expo(i,j)
d(j,2) = ao_extra_coef_normalized(i,j)
enddo
call dsort(d(1,1),iorder,ao_extra_prim_num(i))
call dset_order(d(1,2),iorder,ao_extra_prim_num(i))
do j=1,ao_extra_prim_num(i)
ao_extra_expo_ordered(i,j) = d(j,1)
ao_extra_coef_normalized_ordered(i,j) = d(j,2)
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_extra_coef_normalized_ordered_transp, (ao_extra_prim_num_max,ao_extra_num) ]
implicit none
BEGIN_DOC
! Transposed :c:data:`ao_extra_coef_normalized_ordered`
END_DOC
integer :: i,j
do j=1, ao_extra_num
do i=1, ao_extra_prim_num_max
ao_extra_coef_normalized_ordered_transp(i,j) = ao_extra_coef_normalized_ordered(j,i)
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_extra_expo_ordered_transp, (ao_extra_prim_num_max,ao_extra_num) ]
implicit none
BEGIN_DOC
! Transposed :c:data:`ao_extra_expo_ordered`
END_DOC
integer :: i,j
do j=1, ao_extra_num
do i=1, ao_extra_prim_num_max
ao_extra_expo_ordered_transp(i,j) = ao_extra_expo_ordered(j,i)
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ integer, ao_extra_l, (ao_extra_num) ]
&BEGIN_PROVIDER [ integer, ao_extra_l_max ]
&BEGIN_PROVIDER [ character*(128), ao_extra_l_char, (ao_extra_num) ]
implicit none
BEGIN_DOC
! :math:`l` value of the |ao_extra|: :math`a+b+c` in :math:`x^a y^b z^c`
END_DOC
integer :: i
do i=1,ao_extra_num
ao_extra_l(i) = ao_extra_power(i,1) + ao_extra_power(i,2) + ao_extra_power(i,3)
ao_extra_l_char(i) = l_to_character(ao_extra_l(i))
enddo
ao_extra_l_max = maxval(ao_extra_l)
END_PROVIDER
integer function ao_extra_power_index(nx,ny,nz)
implicit none
integer, intent(in) :: nx, ny, nz
BEGIN_DOC
! Unique index given to a triplet of powers:
!
! :math:`\frac{1}{2} (l-n_x) (l-n_x+1) + n_z + 1`
END_DOC
integer :: l
l = nx + ny + nz
ao_extra_power_index = ((l-nx)*(l-nx+1))/2 + nz + 1
end
BEGIN_PROVIDER [ integer, Nucl_N_ao_extras, (extra_nucl_num)]
&BEGIN_PROVIDER [ integer, N_ao_extras_max ]
implicit none
BEGIN_DOC
! Number of |ao_extras| per atom
END_DOC
integer :: i
Nucl_N_ao_extras = 0
do i = 1, ao_extra_num
Nucl_N_ao_extras(ao_extra_nucl(i)) +=1
enddo
N_ao_extras_max = maxval(Nucl_N_ao_extras)
END_PROVIDER
BEGIN_PROVIDER [ integer, Nucl_ao_extras, (extra_nucl_num,N_ao_extras_max)]
implicit none
BEGIN_DOC
! List of |ao_extras| centered on each atom
END_DOC
integer :: i
integer, allocatable :: nucl_tmp(:)
allocate(nucl_tmp(nucl_num))
nucl_tmp = 0
Nucl_ao_extras = 0
do i = 1, ao_extra_num
nucl_tmp(ao_extra_nucl(i))+=1
Nucl_ao_extras(ao_extra_nucl(i),nucl_tmp(ao_extra_nucl(i))) = i
enddo
deallocate(nucl_tmp)
END_PROVIDER
BEGIN_PROVIDER [ integer, Nucl_list_shell_ao_extras, (extra_nucl_num,N_ao_extras_max)]
&BEGIN_PROVIDER [ integer, Nucl_num_shell_ao_extras, (nucl_num)]
implicit none
integer :: i,j,k
BEGIN_DOC
! Index of the shell type |ao_extras| and of the corresponding |ao_extras|
! By convention, for p,d,f and g |ao_extras|, we take the index
! of the |ao_extra| with the the corresponding power in the x axis
END_DOC
do i = 1, extra_nucl_num
Nucl_num_shell_ao_extras(i) = 0
do j = 1, Nucl_N_ao_extras(i)
if (ao_extra_power(Nucl_ao_extras(i,j),1) == ao_extra_l(Nucl_ao_extras(i,j))) then
Nucl_num_shell_ao_extras(i)+=1
Nucl_list_shell_ao_extras(i,Nucl_num_shell_ao_extras(i))=Nucl_ao_extras(i,j)
endif
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ character*(4), ao_extra_l_char_space, (ao_extra_num) ]
implicit none
BEGIN_DOC
! Converts an l value to a string
END_DOC
integer :: i
character*(4) :: give_ao_extra_character_space
do i=1,ao_extra_num
if(ao_extra_l(i)==0)then
! S type ao_extra
give_ao_extra_character_space = 'S '
elseif(ao_extra_l(i) == 1)then
! P type ao_extra
if(ao_extra_power(i,1)==1)then
give_ao_extra_character_space = 'X '
elseif(ao_extra_power(i,2) == 1)then
give_ao_extra_character_space = 'Y '
else
give_ao_extra_character_space = 'Z '
endif
elseif(ao_extra_l(i) == 2)then
! D type ao_extra
if(ao_extra_power(i,1)==2)then
give_ao_extra_character_space = 'XX '
elseif(ao_extra_power(i,2) == 2)then
give_ao_extra_character_space = 'YY '
elseif(ao_extra_power(i,3) == 2)then
give_ao_extra_character_space = 'ZZ '
elseif(ao_extra_power(i,1) == 1 .and. ao_extra_power(i,2) == 1)then
give_ao_extra_character_space = 'XY '
elseif(ao_extra_power(i,1) == 1 .and. ao_extra_power(i,3) == 1)then
give_ao_extra_character_space = 'XZ '
else
give_ao_extra_character_space = 'YZ '
endif
elseif(ao_extra_l(i) == 3)then
! F type ao_extra
if(ao_extra_power(i,1)==3)then
give_ao_extra_character_space = 'XXX '
elseif(ao_extra_power(i,2) == 3)then
give_ao_extra_character_space = 'YYY '
elseif(ao_extra_power(i,3) == 3)then
give_ao_extra_character_space = 'ZZZ '
elseif(ao_extra_power(i,1) == 2 .and. ao_extra_power(i,2) == 1)then
give_ao_extra_character_space = 'XXY '
elseif(ao_extra_power(i,1) == 2 .and. ao_extra_power(i,3) == 1)then
give_ao_extra_character_space = 'XXZ '
elseif(ao_extra_power(i,2) == 2 .and. ao_extra_power(i,1) == 1)then
give_ao_extra_character_space = 'YYX '
elseif(ao_extra_power(i,2) == 2 .and. ao_extra_power(i,3) == 1)then
give_ao_extra_character_space = 'YYZ '
elseif(ao_extra_power(i,3) == 2 .and. ao_extra_power(i,1) == 1)then
give_ao_extra_character_space = 'ZZX '
elseif(ao_extra_power(i,3) == 2 .and. ao_extra_power(i,2) == 1)then
give_ao_extra_character_space = 'ZZY '
elseif(ao_extra_power(i,3) == 1 .and. ao_extra_power(i,2) == 1 .and. ao_extra_power(i,3) == 1)then
give_ao_extra_character_space = 'XYZ '
endif
elseif(ao_extra_l(i) == 4)then
! G type ao_extra
if(ao_extra_power(i,1)==4)then
give_ao_extra_character_space = 'XXXX'
elseif(ao_extra_power(i,2) == 4)then
give_ao_extra_character_space = 'YYYY'
elseif(ao_extra_power(i,3) == 4)then
give_ao_extra_character_space = 'ZZZZ'
elseif(ao_extra_power(i,1) == 3 .and. ao_extra_power(i,2) == 1)then
give_ao_extra_character_space = 'XXXY'
elseif(ao_extra_power(i,1) == 3 .and. ao_extra_power(i,3) == 1)then
give_ao_extra_character_space = 'XXXZ'
elseif(ao_extra_power(i,2) == 3 .and. ao_extra_power(i,1) == 1)then
give_ao_extra_character_space = 'YYYX'
elseif(ao_extra_power(i,2) == 3 .and. ao_extra_power(i,3) == 1)then
give_ao_extra_character_space = 'YYYZ'
elseif(ao_extra_power(i,3) == 3 .and. ao_extra_power(i,1) == 1)then
give_ao_extra_character_space = 'ZZZX'
elseif(ao_extra_power(i,3) == 3 .and. ao_extra_power(i,2) == 1)then
give_ao_extra_character_space = 'ZZZY'
elseif(ao_extra_power(i,1) == 2 .and. ao_extra_power(i,2) == 2)then
give_ao_extra_character_space = 'XXYY'
elseif(ao_extra_power(i,2) == 2 .and. ao_extra_power(i,3) == 2)then
give_ao_extra_character_space = 'YYZZ'
elseif(ao_extra_power(i,1) == 2 .and. ao_extra_power(i,2) == 1 .and. ao_extra_power(i,3) == 1)then
give_ao_extra_character_space = 'XXYZ'
elseif(ao_extra_power(i,2) == 2 .and. ao_extra_power(i,1) == 1 .and. ao_extra_power(i,3) == 1)then
give_ao_extra_character_space = 'YYXZ'
elseif(ao_extra_power(i,3) == 2 .and. ao_extra_power(i,1) == 1 .and. ao_extra_power(i,2) == 1)then
give_ao_extra_character_space = 'ZZXY'
endif
endif
ao_extra_l_char_space(i) = give_ao_extra_character_space
enddo
END_PROVIDER

View File

@ -0,0 +1,68 @@
! ---
BEGIN_PROVIDER [ integer, Nucl_ao_extras_transposed, (N_ao_extras_max,nucl_num)]
BEGIN_DOC
! List of ao_extras attached on each atom
END_DOC
implicit none
integer :: i
integer, allocatable :: nucl_tmp(:)
allocate(nucl_tmp(nucl_num))
nucl_tmp = 0
do i = 1, ao_extra_num
nucl_tmp(ao_extra_nucl(i)) += 1
Nucl_ao_extras_transposed(nucl_tmp(ao_extra_nucl(i)),ao_extra_nucl(i)) = i
enddo
deallocate(nucl_tmp)
END_PROVIDER
! ---
BEGIN_PROVIDER [double precision, ao_extra_expo_ordered_transp_per_nucl, (ao_extra_prim_num_max,N_ao_extras_max,nucl_num) ]
implicit none
integer :: i,j,k,l
do i = 1, nucl_num
do j = 1,Nucl_N_ao_extras(i)
k = Nucl_ao_extras_transposed(j,i)
do l = 1, ao_extra_prim_num(k)
ao_extra_expo_ordered_transp_per_nucl(l,j,i) = ao_extra_expo_ordered_transp(l,k)
enddo
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ integer, ao_extra_power_ordered_transp_per_nucl, (3,N_ao_extras_max,nucl_num) ]
implicit none
integer :: i,j,k,l
do i = 1, nucl_num
do j = 1,Nucl_N_ao_extras(i)
k = Nucl_ao_extras_transposed(j,i)
do l = 1, 3
ao_extra_power_ordered_transp_per_nucl(l,j,i) = ao_extra_power(k,l)
enddo
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_extra_coef_normalized_ordered_transp_per_nucl, (ao_extra_prim_num_max,N_ao_extras_max,nucl_num) ]
implicit none
integer :: i,j,k,l
do i = 1, nucl_num
do j = 1,Nucl_N_ao_extras(i)
k = Nucl_ao_extras_transposed(j,i)
do l = 1, ao_extra_prim_num(k)
ao_extra_coef_normalized_ordered_transp_per_nucl(l,j,i) = ao_extra_coef_normalized_ordered_transp(l,k)
enddo
enddo
enddo
END_PROVIDER

View File

@ -0,0 +1,87 @@
BEGIN_PROVIDER [ double precision, effective_ao_extra_dm, (ao_extra_num, ao_extra_num)]
implicit none
BEGIN_DOC
! effective density matrix : rho_pq x N_p x N_q x E_pq x (pi/gamma_pq)^3/2
!
! where rho_pq is the usual density matrix
!
! N_p and N_q are the normalization factors associated with the two Gaussians
!
! E_pq is the prefactor resulting from the Gaussian product
!
! gamma_pq = gamma_p + gamm_q
END_DOC
integer :: i,j
do i = 1, ao_extra_num
do j = 1, ao_extra_num
effective_ao_extra_dm(j,i) = ao_extra_one_e_dm(j,i) * ao_extra_coef_normalized(j,1) * ao_extra_coef_normalized(i,1) &
* inv_pi_gamma_pq_3_2_ao_extra(j,i) * E_pq_ao_extra(j,i)
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, gamma_pq_ao_extra, (ao_extra_num,ao_extra_num)]
&BEGIN_PROVIDER [ double precision, inv_gamma_pq_ao_extra, (ao_extra_num,ao_extra_num)]
&BEGIN_PROVIDER [ double precision, inv_pi_gamma_pq_3_2_ao_extra, (ao_extra_num,ao_extra_num)]
&BEGIN_PROVIDER [ double precision, sqrt_gamma_pq_ao_extra, (ao_extra_num,ao_extra_num)]
implicit none
BEGIN_DOC
! gamma_pq_ao_extra = gamma_p + gamma_q
!
! inv_gamma_pq_ao_extra = 1/(gamma_p + gamma_q)
!
! sqrt_gamma_pq_ao_extra = sqrt(gamma_p + gamma_q)
!
! WARNING :: VALID ONLY IN THE CASE OF A PURELY 1S BASIS
END_DOC
include 'constants.include.F'
integer :: i,j
do i = 1, ao_extra_num
do j = 1, ao_extra_num
gamma_pq_ao_extra(j,i) = ao_extra_expo(j,1) + ao_extra_expo(i,1)
inv_gamma_pq_ao_extra(j,i) = 1.d0/gamma_pq_ao_extra(j,i)
sqrt_gamma_pq_ao_extra(j,i) = dsqrt(gamma_pq_ao_extra(j,i))
inv_pi_gamma_pq_3_2_ao_extra(j,i) = (pi * inv_gamma_pq_ao_extra(j,i))**(1.5d0)
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, E_pq_ao_extra, (ao_extra_num,ao_extra_num)]
implicit none
BEGIN_DOC
! E_pq_ao_extra = exp(-alpha_p alpha_q/gamma_pq (Q_p - Q_q)^2)
END_DOC
integer :: i,j
do i = 1, ao_extra_num
do j = 1, ao_extra_num
E_pq_ao_extra(j,i) = dexp(-ao_extra_center_1s_dist(j,i)**2 * ao_extra_expo(j,1)*ao_extra_expo(i,1)*inv_gamma_pq_ao_extra(j,i))
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_extra_center_1s, (3,ao_extra_num)]
implicit none
BEGIN_DOC
! Original position of each extra AO
END_DOC
integer :: i,i_nucl
do i = 1, ao_extra_num
i_nucl= ao_extra_nucl(i)
ao_extra_center_1s(1:3,i) = extra_nucl_coord(i_nucl,1:3)
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_extra_center_1s_dist, (ao_extra_num,ao_extra_num)]
implicit none
integer :: i,j
do i = 1, ao_extra_num
do j = 1, ao_extra_num
ao_extra_center_1s_dist(j,i) = (ao_extra_center_1s(1,j) - ao_extra_center_1s(1,i))**2 &
+ (ao_extra_center_1s(2,j) - ao_extra_center_1s(2,i))**2 &
+ (ao_extra_center_1s(3,j) - ao_extra_center_1s(3,i))**2
ao_extra_center_1s_dist(j,i)=dsqrt(ao_extra_center_1s_dist(j,i))
enddo
enddo
END_PROVIDER

View File

@ -0,0 +1,19 @@
BEGIN_PROVIDER [ integer, n_pt_max_extra_basis_integrals ]
&BEGIN_PROVIDER [ integer, n_pt_max_extra_basis_i_x]
implicit none
BEGIN_DOC
! Number of points used in the numerical integrations.
END_DOC
integer :: n_pt_sup
integer :: prim_power_l_max
include 'utils/constants.include.F'
prim_power_l_max = maxval(ao_extra_power)
n_pt_max_extra_basis_integrals = 24 * prim_power_l_max + 4
n_pt_max_extra_basis_i_x = 8 * prim_power_l_max
ASSERT (n_pt_max_extra_basis_i_x-1 <= max_dim)
if (n_pt_max_extra_basis_i_x-1 > max_dim) then
print *, 'Increase max_dim in utils/constants.include.F to ', n_pt_max_extra_basis_i_x-1
stop 1
endif
END_PROVIDER

View File

@ -0,0 +1,16 @@
program extra_basis
implicit none
BEGIN_DOC
! TODO : Put the documentation of the program here
END_DOC
integer :: i
print*,'extra_nucl_num = ',extra_nucl_num
do i = 1, extra_nucl_num
print*,'i = ',i
print*,'extra_nucl_label = ',extra_nucl_label(i)
print*,'extra_nucl_charge = ',extra_nucl_charge(i)
print*,'extra_nucl_coord = '
print*,extra_nucl_coord(i,1:3)
enddo
print*,ao_extra_num
end

View File

@ -0,0 +1,41 @@
program fit_1s_basis
implicit none
provide lmax_too_big
integer :: i,j
print*,'////////////////////////////////////////////////////'
print*,'////////////////////////////////////////////////////'
print*,'Fitting the original basis set on uncontracted s only functions '
print*,'WARNING :: works for now with only P functions at most !!'
print*,'WARNING :: otherwise it will stop '
print*,'Writting the results in the extra_nuclei and ao_extra_basis folders of EZFIO'
print*,'New number of atomic functions : '
print*,'n_func_tot = ',n_func_tot
print*,'extra_fictious_nucl = ',extra_fictious_nucl
do i = 1, extra_fictious_nucl
print*,list_fict_nucl(i)
enddo
print*,''
print*,''
do i = 1, nucl_num
print*,list_real_nucl(i)
enddo
call ezfio_set_extra_nuclei_extra_nucl_num(new_nucl_num)
call ezfio_set_extra_nuclei_extra_nucl_fictious_num(extra_fictious_nucl)
call ezfio_set_extra_nuclei_extra_nucl_real_num(nucl_num)
call ezfio_set_extra_nuclei_extra_nucl_fictious_list(list_fict_nucl)
call ezfio_set_extra_nuclei_extra_nucl_real_list(list_real_nucl)
call ezfio_set_extra_nuclei_extra_nucl_real_fictious_list(extra_nucl_real_fictious_list_prov)
call ezfio_set_extra_nuclei_extra_nucl_charge(new_nucl_charge_1s)
call ezfio_set_extra_nuclei_extra_nucl_coord(new_nucl_coord_1s)
call ezfio_set_extra_nuclei_extra_nucl_label(new_nucl_label_1s)
!
call ezfio_set_ao_extra_basis_ao_extra_num(n_func_tot)
call ezfio_set_ao_extra_basis_ao_extra_center(ao_extra_center)
call ezfio_set_ao_extra_basis_ao_extra_nucl(new_ao_nucl_1s)
call ezfio_set_ao_extra_basis_ao_extra_prim_num(new_ao_prim_num_1s)
call ezfio_set_ao_extra_basis_ao_extra_coef(new_ao_coef_1s)
call ezfio_set_ao_extra_basis_ao_extra_expo(new_ao_expo_1s)
call ezfio_set_ao_extra_basis_ao_extra_power(new_ao_power_1s)
end

View File

@ -0,0 +1,7 @@
3
O 0.000000 -0.399441 3.000000
H 0.761232 0.199721 3.000000
H -0.761232 0.199721 3.000000

23
src/ao_extra_basis/install Executable file
View File

@ -0,0 +1,23 @@
#!/bin/bash
# Check if the QP_ROOT environment variable is set.
if [[ -z ${QP_ROOT} ]]
then
print "The QP_ROOT environment variable is not set."
print "Please reload the quantum_package.rc file."
exit -1
fi
# Get the absolute path of the current directory.
currdir=${PWD}
# list of the scripts to be used by the module
scripts_list="qp_copy_extra_basis qp_add_extra_fit_system qp_copy_extra_basis_to_usual_basis qp_fit_1s_basis"
# Make a symbolic link for all scripts to be used in the ${QP_ROOT}/scripts/
# directory.
for i in $scripts_list
do
ln --symbolic ${currdir}/$i ${QP_ROOT}/scripts/
done

View File

@ -0,0 +1,266 @@
BEGIN_PROVIDER [ double precision, ao_extra_center]
implicit none
ao_extra_center = 0.01d0
END_PROVIDER
BEGIN_PROVIDER [ integer, n_func_tot]
implicit none
BEGIN_DOC
! n_func_tot :: total number of functions in the fitted basis set
!
! returned in an uncontracted way
END_DOC
integer :: i,prefact
n_func_tot = 0
print*,'n_func_tot '
do i = 1, ao_num
if(ao_l(i) == 0)then
prefact = 1 ! s functions
else
! p functions are fitted with 2 functions
! d functions are fitted with 4 functions etc ...
prefact=2*ao_l(i)
endif
n_func_tot += prefact * ao_prim_num(i)
enddo
END_PROVIDER
BEGIN_PROVIDER [ integer, n_prim_tot_orig]
implicit none
integer :: i
n_prim_tot_orig = 0
do i = 1, ao_num
n_prim_tot_orig += ao_prim_num(i)
enddo
END_PROVIDER
BEGIN_PROVIDER [ logical, lmax_too_big]
implicit none
if (ao_l_max.gt.1)then
lmax_too_big = .True.
else
lmax_too_big = .False.
endif
if(lmax_too_big)then
print*,'STOPPING !! lmax is larger than 1 !'
print*,'Cannot yet fit with 1s functions ...'
stop
endif
END_PROVIDER
BEGIN_PROVIDER [ integer, n_2p_func_orig]
&BEGIN_PROVIDER [ integer, n_2p_func_tot]
implicit none
integer :: i
BEGIN_DOC
! n_2p_func_orig :: number of 2p functions in the original basis
!
! n_2p_func_tot :: total number of p functions in the fitted basis
END_DOC
n_2p_func_orig= 0
n_2p_func_tot = 0
do i = 1, ao_num
if(ao_l(i)==1)then
n_2p_func_orig+= 1
n_2p_func_tot += ao_prim_num(i) * 2
endif
enddo
print*,'n_2p_func_tot = ',n_2p_func_tot
END_PROVIDER
BEGIN_PROVIDER [ integer, list_2p_functions, (n_2p_func_orig)]
implicit none
BEGIN_DOC
! list of 2p functions in the original basis
END_DOC
integer :: i,j
j=0
do i = 1, ao_num
if(ao_l(i)==1)then
j+=1
list_2p_functions(j) = i
endif
enddo
END_PROVIDER
BEGIN_PROVIDER [ integer, extra_fictious_nucl]
implicit none
extra_fictious_nucl = n_2p_func_tot
END_PROVIDER
BEGIN_PROVIDER [ integer, new_nucl_num]
implicit none
new_nucl_num = nucl_num + n_2p_func_tot
print*,'new_nucl_num = ',new_nucl_num
END_PROVIDER
BEGIN_PROVIDER [ character*(32), new_nucl_label_1s , (new_nucl_num) ]
&BEGIN_PROVIDER [ integer, list_real_nucl, (nucl_num) ]
&BEGIN_PROVIDER [ integer, list_fict_nucl, (extra_fictious_nucl) ]
implicit none
integer :: i,j
do i = 1, nucl_num
new_nucl_label_1s(i) = nucl_label(i)
list_real_nucl(i) = i
enddo
j=0
do i = nucl_num+1,new_nucl_num
j+=1
new_nucl_label_1s(i) = "X"
list_fict_nucl(j) = i
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, new_nucl_coord_1s, (new_nucl_num,3)]
implicit none
integer :: i,j
do i = 1, new_nucl_num
new_nucl_coord_1s(i,1:3) = new_nucl_coord_1s_transp(1:3,i)
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, new_nucl_coord_1s_transp, (3,new_nucl_num)]
&BEGIN_PROVIDER [ double precision, new_nucl_charge_1s, (new_nucl_num)]
&BEGIN_PROVIDER [ integer, extra_nucl_real_fictious_list_prov, (extra_fictious_nucl)]
implicit none
BEGIN_DOC
! the real atoms are located in the first nucl_num entries
!
! then the fictious atoms are located after
END_DOC
integer :: i,ii,j,i_ao,k,n_ao
integer :: return_xyz_int,power(3),good_i
new_nucl_charge_1s = 0.d0
do i = 1, nucl_num
new_nucl_coord_1s_transp(1:3,i) = nucl_coord_transp(1:3,i)
new_nucl_charge_1s(i) = nucl_charge(i)
enddo
k = nucl_num
do i = 1, nucl_num
do ii = 1, Nucl_N_Aos(i)
i_ao = nucl_aos_transposed(ii,i)
if(ao_l(i_ao)==1)then
! split the function into 2 s functions
! one is centered in R_x + d
power(1:3) = ao_power(i_ao,1:3)
good_i = return_xyz_int(power)
do j = 1, ao_prim_num(i_ao)
k+=1
new_nucl_coord_1s_transp(1:3,k)= nucl_coord_transp(1:3,i)
new_nucl_coord_1s_transp(good_i,k)+= ao_extra_center
new_nucl_charge_1s(k) = 0.d0
extra_nucl_real_fictious_list_prov(k-nucl_num)=i
k+=1
! one is centered in R_x - d
new_nucl_coord_1s_transp(1:3,k)= nucl_coord_transp(1:3,i)
new_nucl_coord_1s_transp(good_i,k)-= ao_extra_center
new_nucl_charge_1s(k) = 0.d0
extra_nucl_real_fictious_list_prov(k-nucl_num)=i
enddo
else if(ao_l(i_ao).gt.1)then
print*,'WARNING ! Lmax value not implemented yet !'
print*,'stopping ...'
stop
endif
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ integer, new_n_AOs_max]
implicit none
new_n_AOs_max = ao_prim_num_max * n_AOs_max
END_PROVIDER
BEGIN_PROVIDER [ integer, new_Nucl_N_Aos, (new_nucl_num)]
&BEGIN_PROVIDER [ integer, new_nucl_aos_transposed, (new_n_AOs_max,new_nucl_num) ]
&BEGIN_PROVIDER [ double precision, new_ao_expo_1s , (n_func_tot) ]
&BEGIN_PROVIDER [ integer, new_ao_nucl_1s, (n_func_tot)]
implicit none
integer :: i,j,ii,i_ao,n_func,n_func_total,n_nucl
double precision :: coef
n_func_total = 0
do i = 1, nucl_num
n_func = 0
do ii = 1, Nucl_N_Aos(i)
i_ao = nucl_aos_transposed(ii,i)
if(ao_l(i_ao)==0)then
do j = 1, ao_prim_num(i_ao)
coef= ao_expo(i_ao,j)
n_func_total += 1
n_func +=1
new_nucl_aos_transposed(n_func,i) = n_func_total
new_ao_expo_1s(n_func_total) = coef
new_ao_nucl_1s(n_func_total) = i
enddo
endif
enddo
new_Nucl_N_Aos(i) = n_func
enddo
n_nucl=nucl_num
do i = 1, nucl_num
do ii = 1, Nucl_N_Aos(i)
i_ao = nucl_aos_transposed(ii,i)
if(ao_l(i_ao)==1)then
do j = 1, ao_prim_num(i_ao)
coef= ao_expo(i_ao,j)
n_func_total+=1
n_nucl +=1
new_nucl_aos_transposed(1,n_nucl) = n_func_total
new_ao_expo_1s(n_func_total) = coef
new_Nucl_N_Aos(n_nucl)=1
new_ao_nucl_1s(n_func_total) = n_nucl
n_func_total+=1
n_nucl +=1
new_nucl_aos_transposed(1,n_nucl) = n_func_total
new_ao_expo_1s(n_func_total) = coef
new_Nucl_N_Aos(n_nucl)=1
new_ao_nucl_1s(n_func_total) = n_nucl
enddo
endif
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, new_ao_coef_1s , (n_func_tot) ]
implicit none
integer :: i
BEGIN_DOC
! Primitive coefficients, read from input. Those should not be used directly, as the MOs are expressed on the basis of **normalized** AOs.
END_DOC
do i = 1, n_func_tot
new_ao_coef_1s(i) = 1.d0
enddo
END_PROVIDER
BEGIN_PROVIDER [ integer, new_ao_prim_num_1s, (n_func_tot)]
implicit none
integer :: i
do i = 1, n_func_tot
new_ao_prim_num_1s(i) = 1
enddo
END_PROVIDER
BEGIN_PROVIDER [integer, new_ao_power_1s, (n_func_tot,3)]
implicit none
new_ao_power_1s = 0
END_PROVIDER
integer function return_xyz_int(power)
implicit none
integer,intent(in) :: power(3)
if(power(1) == 1 .and. power(2) ==0 .and. power(3) ==0)then
return_xyz_int = 1
else if (power(2) == 1 .and. power(1) ==0 .and. power(3) ==0)then
return_xyz_int = 2
else if (power(3) == 1 .and. power(1) ==0 .and. power(1) ==0)then
return_xyz_int = 3
else
return_xyz_int = -1000
endif
end

View File

@ -0,0 +1,19 @@
#!/bin/bash
# specify the QP folder
QP=$QP_ROOT
dir=${QP}
# sourcing the quantum_package.rc file
. ${QP}/quantum_package.rc
# The main EZFIO folder
main=${1%/}
# The XYZ file containing the geometry of the additional system you want to add
extra=${2%.xyz}
basis_extra=sto-3g
ezfio_extra=${extra}_${basis_extra}_1s
echo $ezfio_extra
qp_fit_1s_basis $extra $basis_extra
qp set_file $ezfio_extra
qp run scf | tee ${ezfio_extra}.scf.out
qp run save_one_e_dm | tee ${ezfio_extra}.one_rdm.out
qp_copy_extra_basis ${ezfio_extra} $main

View File

@ -0,0 +1,67 @@
#!/bin/bash
# specify the QP folder
QP=$QP_ROOT
dir=${QP}
# sourcing the quantum_package.rc file
. ${QP}/quantum_package.rc
# script that copy all data from |AO| basis and nuclei of EZFIO_extra to the ao_extra_basis and extra_nuclei of EZFIO_target
# use:
# qp_copy_extra_basis EZFIO_extra EZFIO_target
EZFIO_extra=${1%/}
EZFIO_extra=${EZFIO_extra%.xyz}
EZFIO_target=${2%/}
echo "********** SCRIPT TO COPY DATA FROM EZFIO TO ANOTHER *********"
echo "Extracting data from "$EZFIO_extra
echo "and Copying data to "$EZFIO_target
### COPYING ALL DATA FROM $EZFIO_extra/nuclei/ to ${EZFIO_target}/extra_nuclei/
echo "COPYING ALL DATA FROM "$EZFIO_extra"/nuclei/ to "${EZFIO_target}"/extra_nuclei/"
direxists=false
if [ -d ${EZFIO_target}/extra_nuclei/ ] ; then
direxists=true
echo "The directory extra_nuclei exists"
else
echo "Creating the directory extra_nuclei "
direxists=false
mkdir ${EZFIO_target}/extra_nuclei/
fi
data=`\ls $EZFIO_extra/nuclei/`
for i in $data
do
echo $i
newfile=`echo $i | sed 's/nucl/extra_nucl/g' `
echo $newfile
cp ${EZFIO_extra}/nuclei/$i ${EZFIO_target}/extra_nuclei/$newfile
done
### COPYING ALL DATA FROM $EZFIO_extra/ao_basis/ to ${EZFIO_target}/ao_extra_basis/
direxists=false
if [ -d ${EZFIO_target}/ao_extra_basis/ ] ; then
direxists=true
echo "The directory exists ao_extra_basis"
else
echo "Creating the directory ao_extra_basis"
direxists=false
mkdir ${EZFIO_target}/ao_extra_basis/
fi
echo "COPYING ALL DATA FROM "$EZFIO_extra"/ao_basis/ to "${EZFIO_target}"/ao_extra_basis/"
data=`\ls $EZFIO_extra/ao_basis/`
for i in $data
do
echo $i
newfile=`echo $i | sed 's/ao/ao_extra/g' `
echo $newfile
cp ${EZFIO_extra}/ao_basis/$i ${EZFIO_target}/ao_extra_basis/$newfile
done
i=primitives_normalized
newfile=primitives_normalized_extra
cp ${EZFIO_extra}/ao_basis/$i ${EZFIO_target}/ao_extra_basis/$newfile
echo "COPYING ALL DATA FROM "$EZFIO_extra"/aux_quantities/ to "${EZFIO_target}"/ao_extra_basis/"
i=data_one_e_dm_tot_ao.gz
newfile=ao_extra_one_e_dm.gz
cp ${EZFIO_extra}/aux_quantities/$i ${EZFIO_target}/ao_extra_basis/$newfile

View File

@ -0,0 +1,57 @@
#!/bin/bash
# specify the QP folder
# script that copy all data in ao_extra_basis and extra_nuclei and copy it to ao_basis and ao_nuclei
QP=$QP_ROOT
dir=${QP}
# sourcing the quantum_package.rc file
. ${QP}/quantum_package.rc
EZFIO_extra=${1%/}
EZFIO_target=${2%/}
echo "********** SCRIPT TO COPY DATA FROM EZFIO TO ANOTHER *********"
echo "Extracting data from "$EZFIO_extra
echo "and Copying data to "$EZFIO_target
### COPYING ALL DATA FROM $EZFIO_extra/nuclei/ to ${EZFIO_target}/extra_nuclei/
echo "COPYING ALL DATA FROM "$EZFIO_extra"/extra_nuclei/ to "${EZFIO_target}"/nuclei/"
direxists=false
if [ -d ${EZFIO_target}/extra_nuclei/ ] ; then
direxists=true
echo "The directory extra_nuclei exists"
else
echo "PB !!"
direxists=false
mkdir ${EZFIO_target}/extra_nuclei/
fi
data=`\ls $EZFIO_extra/extra_nuclei/`
for i in $data
do
echo $i
newfile=`echo $i | sed 's/extra_nucl/nucl/g' `
echo $newfile
cp ${EZFIO_extra}/extra_nuclei/$i ${EZFIO_target}/nuclei/$newfile
done
### COPYING ALL DATA FROM $EZFIO_extra/ao_basis/ to ${EZFIO_target}/ao_extra_basis/
direxists=false
if [ -d ${EZFIO_target}/ao_extra_basis/ ] ; then
direxists=true
echo "The directory exists ao_extra_basis"
else
echo "PB !!"
direxists=false
mkdir ${EZFIO_target}/ao_extra_basis/
fi
echo "COPYING ALL DATA FROM "$EZFIO_extra"/ao_basis/ to "${EZFIO_target}"/ao_extra_basis/"
data=`\ls $EZFIO_extra/ao_extra_basis/`
for i in $data
do
echo $i
newfile=`echo $i | sed 's/ao_extra/ao/g' `
echo $newfile
cp ${EZFIO_extra}/ao_extra_basis/$i ${EZFIO_target}/ao_basis/$newfile
done

View File

@ -0,0 +1,12 @@
#!/bin/bash
## Takes as an argument an xyz file and a basis, and fit the AO basis onto an "s" basis only
## use:
# qp_fit_1s_basis my_xyz_file.xyz basis
source ~/qp2/quantum_package.rc
input=${1%.xyz}
basis=$2
ezfio_fit=${input}_${basis}_1s
qp create_ezfio -b $basis $input.xyz -o $ezfio_fit
# Fitting the original basis on 1s only basis functions
qp run fit_1s_basis |tee ${ezfio_fit}.fit_1s.out
qp_copy_extra_basis_to_usual_basis $ezfio_fit $ezfio_fit

View File

@ -0,0 +1,49 @@
subroutine rotate_nuclei(phi,theta,psi,nucl_centers,n_nucl,nucl_centers_after)
implicit none
BEGIN_DOC
! routine that rotates a set of nuclei according to three axis corresponding to angles phi, theta, psi.
END_DOC
double precision, intent(in) :: phi,theta,psi
double precision, intent(in) :: nucl_centers(3,n_nucl)
integer, intent(in) :: n_nucl
double precision, intent(out):: nucl_centers_after(3,n_nucl)
double precision :: r_mat(3,3)
call r_phi_theta_psi_matrix(phi,theta,psi,r_mat)
call get_AB_prod(r_mat,3,3,nucl_centers,n_nucl,nucl_centers_after)
end
subroutine r_phi_theta_psi_matrix(phi,theta,psi,r_mat)
implicit none
BEGIN_DOC
! routine that creates the rotation matrix corresponding to phi,theta,psi
!
! according to conventions in MDFT code
END_DOC
double precision, intent(in) :: phi,theta,psi
double precision, intent(out):: r_mat(3,3)
double precision :: ctheta, stheta
double precision :: cphi , sphi
double precision :: cpsi , spsi
ctheta = dcos(theta)
cphi = dcos(phi)
cpsi = dcos(psi)
stheta = dsin(theta)
sphi = dsin(phi)
spsi = dsin(psi)
r_mat(1,1) = ctheta*cphi*cpsi-sphi*spsi
r_mat(1,2) = -ctheta*cphi*spsi-sphi*cpsi
r_mat(1,3) = stheta*cphi
r_mat(2,1) = ctheta*sphi*cpsi+cphi*spsi
r_mat(2,2) = -ctheta*sphi*spsi+cphi*cpsi
r_mat(2,3) = stheta*sphi
r_mat(3,1) = -stheta*cpsi
r_mat(3,2) = stheta*spsi
r_mat(3,3) = ctheta
end

20
src/ao_extra_basis/uninstall Executable file
View File

@ -0,0 +1,20 @@
#!/bin/bash
# Check if the QP_ROOT environment variable is set.
if [[ -z ${QP_ROOT} ]]
then
print "The QP_ROOT environment variable is not set."
print "Please reload the quantum_package.rc file."
exit -1
fi
# list of the scripts to be used by the module
scripts_list="qp_copy_extra_basis qp_add_extra_fit_system qp_copy_extra_basis_to_usual_basis qp_fit_1s_basis"
# Destroy ONLY the symbolic link for the scripts to be used in the
# ${QP_ROOT}/scripts/ directory.
for i in $scripts_list
do
find ${QP_ROOT}/scripts/$i -type l -delete
done

View File

@ -0,0 +1,56 @@
[extra_nucl_num]
doc: Number of nuclei
type: integer
interface: ezfio, provider
[extra_nucl_pouet]
doc: Number of nuclei
type: integer
interface: ezfio, provider, ocaml
default:1
[extra_nucl_label]
doc: Nuclear labels
type: character*(32)
size: (extra_nuclei.extra_nucl_num)
interface: ezfio, provider
[extra_nucl_charge]
doc: Nuclear charges
type:double precision
size: (extra_nuclei.extra_nucl_num)
interface: ezfio, provider
[extra_nucl_coord]
doc: Nuclear coordinates in the format (:, {x,y,z})
type: double precision
size: (extra_nuclei.extra_nucl_num,3)
interface: ezfio
[extra_nucl_real_num]
doc: Number of real nuclei
type: integer
interface: ezfio, provider
[extra_nucl_fictious_num]
doc: Number of fictious nuclei
type: integer
interface: ezfio, provider
[extra_nucl_real_fictious_list]
doc: List of real nuclei to which fictious nuclei are attached to
type: integer
size: (extra_nuclei.extra_nucl_fictious_num)
interface: ezfio, provider
[extra_nucl_fictious_list]
doc: List of fictious nuclei
type: integer
size: (extra_nuclei.extra_nucl_fictious_num)
interface: ezfio, provider
[extra_nucl_real_list]
doc: List of real nuclei
type: integer
size: (extra_nuclei.extra_nucl_real_num)
interface: ezfio, provider

3
src/extra_nuclei/NEED Normal file
View File

@ -0,0 +1,3 @@
ezfio_files
utils
nuclei

View File

@ -0,0 +1,4 @@
============
extra_nuclei
============

View File

@ -0,0 +1,7 @@
program extra_nuclei
implicit none
BEGIN_DOC
! TODO : Put the documentation of the program here
END_DOC
print *, 'Hello world'
end

View File

@ -0,0 +1,138 @@
BEGIN_PROVIDER [ double precision, extra_nucl_coord, (extra_nucl_num,3) ]
implicit none
BEGIN_DOC
! Nuclear coordinates in the format (:, {x,y,z})
END_DOC
PROVIDE ezfio_filename extra_nucl_label extra_nucl_charge
if (mpi_master) then
double precision, allocatable :: buffer(:,:)
extra_nucl_coord = 0.d0
allocate (buffer(extra_nucl_num,3))
buffer = 0.d0
logical :: has
call ezfio_has_extra_nuclei_extra_nucl_coord(has)
if (.not.has) then
print *, irp_here
stop 1
endif
call ezfio_get_extra_nuclei_extra_nucl_coord(buffer)
integer :: i,j
do i=1,3
do j=1,extra_nucl_num
extra_nucl_coord(j,i) = buffer(j,i)
enddo
enddo
deallocate(buffer)
character*(64), parameter :: f = '(A16, 4(1X,F12.6))'
character*(64), parameter :: ft= '(A16, 4(1X,A12 ))'
double precision, parameter :: a0= 0.529177249d0
call write_time(6)
write(6,'(A)') ''
write(6,'(A)') 'Extra Nuclear Coordinates (Angstroms)'
write(6,'(A)') '====================================='
write(6,'(A)') ''
write(6,ft) &
'================','============','============','============','============'
write(6,*) &
' Atom Charge X Y Z '
write(6,ft) &
'================','============','============','============','============'
do i=1,extra_nucl_num
write(6,f) extra_nucl_label(i), extra_nucl_charge(i), &
extra_nucl_coord(i,1)*a0, &
extra_nucl_coord(i,2)*a0, &
extra_nucl_coord(i,3)*a0
enddo
write(6,ft) &
'================','============','============','============','============'
write(6,'(A)') ''
if (extra_nucl_num > 1) then
double precision :: dist_min, x, y, z
dist_min = huge(1.d0)
do i=1,extra_nucl_num
do j=i+1,extra_nucl_num
x = extra_nucl_coord(i,1)-extra_nucl_coord(j,1)
y = extra_nucl_coord(i,2)-extra_nucl_coord(j,2)
z = extra_nucl_coord(i,3)-extra_nucl_coord(j,3)
dist_min = min(x*x + y*y + z*z, dist_min)
enddo
enddo
write(6,'(A,F12.4,A)') 'Minimal interatomic distance found: ', &
dsqrt(dist_min)*a0,' Angstrom'
endif
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( extra_nucl_coord, 3*extra_nucl_num, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
if (ierr /= MPI_SUCCESS) then
stop 'Unable to read nucl_coord with MPI'
endif
IRP_ENDIF
END_PROVIDER
BEGIN_PROVIDER [ double precision, extra_nucl_coord_transp, (3,extra_nucl_num) ]
implicit none
BEGIN_DOC
! Transposed array of extra_nucl_coord
END_DOC
integer :: i, k
extra_nucl_coord_transp = 0.d0
do i=1,extra_nucl_num
extra_nucl_coord_transp(1,i) = extra_nucl_coord(i,1)
extra_nucl_coord_transp(2,i) = extra_nucl_coord(i,2)
extra_nucl_coord_transp(3,i) = extra_nucl_coord(i,3)
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, extra_center_of_mass, (3) ]
implicit none
BEGIN_DOC
! Center of mass of the molecule
END_DOC
integer :: i,j
double precision :: s
extra_center_of_mass(:) = 0.d0
s = 0.d0
do i=1,extra_nucl_num
do j=1,3
extra_center_of_mass(j) += extra_nucl_coord(i,j)* element_mass(int(extra_nucl_charge(i)))
enddo
s += element_mass(int(extra_nucl_charge(i)))
enddo
s = 1.d0/s
extra_center_of_mass(:) = extra_center_of_mass(:)*s
END_PROVIDER
BEGIN_PROVIDER [ double precision, extra_nucl_dist, (extra_nucl_num,extra_nucl_num)]
implicit none
integer :: i,j
double precision :: x,y,z
do i = 1, extra_nucl_num
do j = 1, extra_nucl_num
x = extra_nucl_coord(i,1)-extra_nucl_coord(j,1)
y = extra_nucl_coord(i,2)-extra_nucl_coord(j,2)
z = extra_nucl_coord(i,3)-extra_nucl_coord(j,3)
extra_nucl_dist(j,i) = x*x+y*y+z*z
extra_nucl_dist(j,i) = dsqrt(extra_nucl_dist(j,i))
enddo
enddo
END_PROVIDER

View File

@ -22,13 +22,3 @@ BEGIN_PROVIDER [ double precision, mo_one_e_integrals,(mo_num,mo_num)]
END_PROVIDER END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_one_e_integrals_from_mo, (ao_num, ao_num)]
implicit none
BEGIN_DOC
! Integrals of the one e hamiltonian obtained from the integrals on the MO basis
!
! WARNING : this is equal to ao_one_e_integrals only if the AO and MO basis have the same number of functions
END_DOC
call mo_to_ao(mo_one_e_integrals,mo_num,ao_one_e_integrals_from_mo,ao_num)
END_PROVIDER