LCPQ/qmcchem
10
1
Fork 1
mirror of https://gitlab.com/scemama/qmcchem.git synced 2024-11-07 06:33:38 +01:00
Code Releases Activity

Added some .irp.f files

Browse Source
This commit is contained in:
Anthony Scemama 2015-12-19 03:29:52 +01:00
parent 8542c1a110
commit ea185b405f
17 changed files with 4558 additions and 3 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

10
README.md
View File

@ -26,9 +26,13 @@ diffusion Monte Carlo algorithm.
in Cloud environments (rance Grilles) coupled to supercomputers in Cloud environments (rance Grilles) coupled to supercomputers
*Warning*: QMC=Chem is under the GPLv2 license. Any modifications to or Warnings:
software including (via compiler) GPL-licensed code must also be made available * QMC=Chem is under the GPLv2 license. Any modifications to or
under the GPL along with build & install instructions. software including (via compiler) GPL-licensed code must also be made available
under the GPL along with build & install instructions.
* Pseudopotentials are about to change in EZFIO database. Current calculations
will not be compatible with future versions
Requirements Requirements
------------ ------------

11
src/constants.F Normal file
View File

@ -0,0 +1,11 @@
real, parameter :: pi=3.14159265359
real, parameter :: two_pi=3.14159265359*2.
real, parameter :: sqpi=1.77245385091
real, parameter :: sq3=1.7320508075688772
real, parameter :: one_over_sqpi = 0.5641895835477563
double precision, parameter :: dpi=3.14159265359d0
double precision, parameter :: dsqpi=1.77245385091d0
double precision, parameter :: dsq3=1.7320508075688772d0
double precision, parameter :: dtwo_pi=3.14159265359d0*2.d0
double precision, parameter :: dfour_pi=3.14159265359d0*4.d0

1551
src/det.irp.f Normal file
View File

File diff suppressed because it is too large Load Diff

333
src/electrons.irp.f Normal file
View File

@ -0,0 +1,333 @@
BEGIN_PROVIDER [ double precision, xbrown, (elec_num_8,3) ]
BEGIN_DOC
! Brownian step. Built in Brownian_step subroutine.
END_DOC
integer :: i,l
xbrown = 0.d0
END_PROVIDER
BEGIN_PROVIDER [ integer, elec_alpha_num ]
&BEGIN_PROVIDER [ integer, elec_alpha_num_8 ]
implicit none
BEGIN_DOC
! Number of alpha electrons
END_DOC
integer, external :: mod_align
elec_alpha_num = -1
call get_electrons_elec_alpha_num(elec_alpha_num)
if (elec_alpha_num <= 0) then
call abrt(irp_here,'Number of alpha electrons should be > 0')
endif
elec_alpha_num_8 = mod_align(elec_alpha_num)
END_PROVIDER
BEGIN_PROVIDER [ integer, elec_beta_num ]
&BEGIN_PROVIDER [ integer, elec_beta_num_8 ]
implicit none
BEGIN_DOC
! Number of beta electrons
END_DOC
integer, external :: mod_align
elec_beta_num = 0
call get_electrons_elec_beta_num(elec_beta_num)
if (elec_beta_num < 0) then
call abrt(irp_here,'Number of beta electrons should be >= 0')
endif
elec_beta_num_8 = mod_align(elec_beta_num)
END_PROVIDER
BEGIN_PROVIDER [ integer, elec_num ]
&BEGIN_PROVIDER [ integer, elec_num_8 ]
&BEGIN_PROVIDER [ integer, elec_num_1_8 ]
implicit none
BEGIN_DOC
! Number of electrons
END_DOC
integer, external :: mod_align
elec_num = elec_alpha_num + elec_beta_num
ASSERT ( elec_num > 0 )
elec_num_8 = mod_align(elec_num)
elec_num_1_8 = mod_align(elec_num+1)
END_PROVIDER
BEGIN_PROVIDER [ real, elec_coord_full, (elec_num+1,3,walk_num) ]
implicit none
BEGIN_DOC
! Electron coordinates of all walkers
! Component (elec_num+1,1,walk_num) contains the length realized by the walker.
! Initialized in init_walkers
END_DOC
integer :: i,k
real, allocatable :: buffer2(:,:,:)
if ( is_worker ) then
call get_elec_coord_full(elec_coord_full,size(elec_coord_full,1))
else
if (.not.do_prepare) then
allocate ( buffer2(elec_num+1,3,elec_coord_pool_size) )
call get_electrons_elec_coord_pool(buffer2)
do k=1,walk_num
do i=1,elec_num+1
elec_coord_full(i,1,k) = buffer2(i,1,k)
elec_coord_full(i,2,k) = buffer2(i,2,k)
elec_coord_full(i,3,k) = buffer2(i,3,k)
enddo
enddo
deallocate ( buffer2 )
else
elec_coord_full = 0.
endif
endif
END_PROVIDER
BEGIN_PROVIDER [ integer, elec_coord_pool_size ]
implicit none
BEGIN_DOC
! Size of the pool of electron coordinates
END_DOC
elec_coord_pool_size = walk_num_tot
call get_electrons_elec_coord_pool_size(elec_coord_pool_size)
call iinfo(irp_here,'elec_coord_pool_size',elec_coord_pool_size)
END_PROVIDER
BEGIN_PROVIDER [ real, elec_coord, (elec_num_1_8,3) ]
implicit none
BEGIN_DOC
! Electron coordinates
END_DOC
integer :: i, k
elec_coord = 0.
do k=1,3
do i=1,elec_num+1
elec_coord(i,k) = elec_coord_full(i,k,walk_i)
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ real, elec_coord_transp, (8,elec_num)
implicit none
BEGIN_DOC
! Transposed array of elec_coord
END_DOC
integer :: i, k
integer, save :: ifirst = 0
if (ifirst == 0) then
ifirst = 1
elec_coord_transp = 0.
endif
!DIR$ VECTOR ALIGNED
!DIR$ LOOP COUNT (200)
do i=1,elec_num
elec_coord_transp(1,i) = elec_coord(i,1)
elec_coord_transp(2,i) = elec_coord(i,2)
elec_coord_transp(3,i) = elec_coord(i,3)
enddo
END_PROVIDER
BEGIN_PROVIDER [ real, elec_dist, (elec_num_8,elec_num) ]
&BEGIN_PROVIDER [ real, elec_dist_vec_x, (elec_num_8,elec_num) ]
&BEGIN_PROVIDER [ real, elec_dist_vec_y, ((-simd_sp+1):elec_num_8,elec_num) ]
&BEGIN_PROVIDER [ real, elec_dist_vec_z, ((-2*simd_sp+1):elec_num_8,elec_num) ]
implicit none
BEGIN_DOC
! Electron-electron distances
END_DOC
integer :: ie1, ie2, l
integer, save :: ifirst = 0
if (ifirst == 0) then
ifirst = 1
!DIR$ VECTOR ALIGNED
elec_dist = 0.
!DIR$ VECTOR ALIGNED
elec_dist_vec_x = 0.
!DIR$ VECTOR ALIGNED
elec_dist_vec_y = 0.
!DIR$ VECTOR ALIGNED
elec_dist_vec_z = 0.
endif
do ie2 = 1,elec_num
real :: x, y, z
x = elec_coord(ie2,1)
y = elec_coord(ie2,2)
z = elec_coord(ie2,3)
!DIR$ VECTOR ALIGNED
!DIR$ LOOP COUNT (200)
do ie1 = 1,elec_num
elec_dist_vec_x(ie1,ie2) = elec_coord(ie1,1) - x
elec_dist_vec_y(ie1,ie2) = elec_coord(ie1,2) - y
elec_dist_vec_z(ie1,ie2) = elec_coord(ie1,3) - z
enddo
!DIR$ VECTOR ALIGNED
!DIR$ LOOP COUNT (200)
do ie1 = 1,elec_num
elec_dist(ie1,ie2) = sqrt( &
elec_dist_vec_x(ie1,ie2)*elec_dist_vec_x(ie1,ie2) + &
elec_dist_vec_y(ie1,ie2)*elec_dist_vec_y(ie1,ie2) + &
elec_dist_vec_z(ie1,ie2)*elec_dist_vec_z(ie1,ie2) )
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ real, nucl_elec_dist, (nucl_num_8,elec_num) ]
&BEGIN_PROVIDER [ real, nucl_elec_dist_vec, (3,nucl_num,elec_num) ]
implicit none
BEGIN_DOC
! Electron-nucleus distances |r_elec - R_nucl|
END_DOC
integer :: i,j,l
integer, save :: ifirst = 0
if (ifirst == 0) then
ifirst = 1
!DIR$ VECTOR ALIGNED
nucl_elec_dist = 0.
!DIR$ VECTOR ALIGNED
nucl_elec_dist_vec = 0.
endif
do i = 1,elec_num
!DIR$ VECTOR ALIGNED
!DIR$ LOOP COUNT (100)
do j = 1,nucl_num
nucl_elec_dist_vec(1,j,i) = elec_coord_transp(1,i) - nucl_coord(j,1)
nucl_elec_dist_vec(2,j,i) = elec_coord_transp(2,i) - nucl_coord(j,2)
nucl_elec_dist_vec(3,j,i) = elec_coord_transp(3,i) - nucl_coord(j,3)
enddo
enddo
do i = 1,elec_num
!DIR$ VECTOR ALIGNED
!DIR$ LOOP COUNT (100)
do j = 1,nucl_num
nucl_elec_dist(j,i) = (elec_coord(i,1) - nucl_coord(j,1)) &
* (elec_coord(i,1) - nucl_coord(j,1)) &
+ (elec_coord(i,2) - nucl_coord(j,2)) &
* (elec_coord(i,2) - nucl_coord(j,2)) &
+ (elec_coord(i,3) - nucl_coord(j,3)) &
* (elec_coord(i,3) - nucl_coord(j,3))
nucl_elec_dist(j,i) = max(1.e-6,sqrt(nucl_elec_dist(j,i)))
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ integer, elec_num_2, (2) ]
BEGIN_DOC
! Number of alpha and beta electrons in an array
END_DOC
elec_num_2(1) = elec_alpha_num
elec_num_2(2) = elec_beta_num
END_PROVIDER
BEGIN_PROVIDER [ integer, elec_spin, (elec_num) ]
implicit none
BEGIN_DOC
! Electron spin. +1 for alpha and -1 for beta
END_DOC
integer :: i
do i=1,elec_alpha_num
elec_spin(i) = 1
enddo
do i=elec_alpha_num+1,elec_num
elec_spin(i) = -1
enddo
END_PROVIDER
BEGIN_PROVIDER [ real, elec_dist_inv, (elec_num_8,elec_num) ]
implicit none
BEGIN_DOC
! 1/rij matrix
END_DOC
integer :: i,j
integer, save :: ifirst = 0
if (ifirst == 0) then
ifirst = 1
elec_dist_inv = 0.
endif
do i=1,elec_num
!DIR$ VECTOR ALIGNED
!DIR$ LOOP COUNT (200)
do j=1,elec_num
elec_dist_inv(j,i) = 1./(elec_dist(j,i)+1.e-12)
enddo
elec_dist_inv(i,i) = 0.
enddo
END_PROVIDER
BEGIN_PROVIDER [ real, nucl_elec_dist_inv, (nucl_num_8,elec_num) ]
implicit none
BEGIN_DOC
! 1/rij matrix
END_DOC
integer :: i,j
do j=1,elec_num
!DIR$ VECTOR ALIGNED
!DIR$ LOOP COUNT (100)
do i=1,nucl_num
nucl_elec_dist_inv(i,j) = 1./nucl_elec_dist(i,j)
enddo
enddo
END_PROVIDER
subroutine save_elec_coord_full
implicit none
BEGIN_DOC
! Save the electron coordinates to disk
END_DOC
integer :: i,k,l
real, allocatable :: buffer2(:,:,:)
allocate ( buffer2(elec_num+1,3,elec_coord_pool_size) )
k=0
do l=1,elec_coord_pool_size
k = k+1
if (k == walk_num+1) then
k=1
endif
do i=1,elec_num+1
buffer2(i,1,l) = elec_coord_full(i,1,k)
buffer2(i,2,l) = elec_coord_full(i,2,k)
buffer2(i,3,l) = elec_coord_full(i,3,k)
enddo
enddo
call ezfio_set_electrons_elec_coord_pool(buffer2)
deallocate ( buffer2 )
end

182
src/ezfio_interface.irp.f Normal file
View File

@ -0,0 +1,182 @@
BEGIN_SHELL [ /usr/bin/python ]
data = [ \
("electrons_elec_coord_pool_size" , "integer" , "" ),
("electrons_elec_coord_pool" , "real" , "(elec_num+1,3,elec_coord_pool_size)" ),
("nuclei_nucl_num" , "integer" , "" ),
("nuclei_nucl_charge" , "real" , "(nucl_num)" ),
("nuclei_nucl_coord" , "real" , "(nucl_num,3)" ),
("nuclei_nucl_fitcusp_radius" , "real" , "(nucl_num)" ),
("mo_basis_mo_coef" , "real" , "(ao_num,mo_tot_num)" ),
("electrons_elec_fitcusp_radius" , "real" , "" ),
("electrons_elec_alpha_num" , "integer" , "" ),
("electrons_elec_beta_num" , "integer" , "" ),
("electrons_elec_walk_num" , "integer" , "" ),
("electrons_elec_walk_num_tot" , "integer" , "" ),
("ao_basis_ao_num" , "integer" , "" ),
("ao_basis_ao_prim_num" , "integer" , "(ao_num)" ),
("ao_basis_ao_nucl" , "integer" , "(ao_num)" ),
("ao_basis_ao_power" , "integer" , "(ao_num,3)" ),
("ao_basis_ao_expo" , "real" , "(ao_num,ao_prim_num_max)" ),
("ao_basis_ao_coef" , "real" , "(ao_num,ao_prim_num_max)" ),
("jastrow_jast_a_up_up" , "real" , "" ),
("jastrow_jast_a_up_dn" , "real" , "" ),
("jastrow_jast_b_up_up" , "real" , "" ),
("jastrow_jast_b_up_dn" , "real" , "" ),
("jastrow_jast_pen" , "real" , "(nucl_num)" ),
("jastrow_jast_eeN_e_a" , "real" , "" ),
("jastrow_jast_eeN_e_b" , "real" , "" ),
("jastrow_jast_eeN_N" , "real" , "(nucl_num)" ),
("jastrow_jast_core_a1" , "real" , "(nucl_num)" ),
("jastrow_jast_core_a2" , "real" , "(nucl_num)" ),
("jastrow_jast_core_b1" , "real" , "(nucl_num)" ),
("jastrow_jast_core_b2" , "real" , "(nucl_num)" ),
("jastrow_jast_type" , "character*(32)", "" ),
("simulation_stop_time" , "integer" , "" ),
("simulation_equilibration" , "logical" , "" ),
("simulation_block_time" , "integer" , "" ),
("simulation_time_step" , "real" , "" ),
("simulation_method" , "character*(32)", "" ),
("simulation_save_data" , "logical" , "" ),
("simulation_print_level" , "integer" , "" ),
("simulation_do_nucl_fitcusp" , "logical" , "" ),
("simulation_sampling" , "character*(32)", "" ),
("simulation_ci_threshold" , "double precision" , "" ),
("simulation_http_server" , "character*(128)", "" ),
("simulation_md5_key" , "character*(32)" , "" ),
("simulation_e_ref" , "double precision" , "" ),
("simulation_do_run" , "logical " , "" ),
("pseudo_do_pseudo" , "logical " , "" ),
]
data_no_set = [\
("mo_basis_mo_tot_num" , "integer" , ""),
("mo_basis_mo_active_num" , "integer" , ""),
("mo_basis_mo_closed_num" , "integer" , ""),
("pseudo_ao_pseudo_grid" , "double precision" , "(ao_num,pseudo_lmax+pseudo_lmax+1,pseudo_lmax-0+1,nucl_num,pseudo_grid_size)"),
("pseudo_mo_pseudo_grid" , "double precision" , "(ao_num,pseudo_lmax+pseudo_lmax+1,pseudo_lmax-0+1,nucl_num,pseudo_grid_size)"),
("pseudo_pseudo_dz_k" , "double precision" , "(nucl_num,pseudo_klocmax)"),
("pseudo_pseudo_dz_kl" , "double precision" , "(nucl_num,pseudo_kmax,pseudo_lmax+1)"),
("pseudo_pseudo_grid_rmax" , "double precision" , ""),
("pseudo_pseudo_grid_size" , "integer" , ""),
("pseudo_pseudo_klocmax" , "integer" , ""),
("pseudo_pseudo_kmax" , "integer" , ""),
("pseudo_pseudo_lmax" , "integer" , ""),
("pseudo_pseudo_n_k" , "integer" , "(nucl_num,pseudo_klocmax)"),
("pseudo_pseudo_n_kl" , "integer" , "(nucl_num,pseudo_kmax,pseudo_lmax+1)"),
("pseudo_pseudo_v_k" , "double precision" , "(nucl_num,pseudo_klocmax)"),
("pseudo_pseudo_v_kl" , "double precision" , "(nucl_num,pseudo_kmax,pseudo_lmax+1)"),
("spindeterminants_n_det_alpha" , "integer" , ""),
("spindeterminants_n_det_beta" , "integer" , ""),
("spindeterminants_n_det" , "integer" , ""),
("spindeterminants_n_int" , "integer" , ""),
("spindeterminants_bit_kind" , "integer" , ""),
("spindeterminants_n_states" , "integer" , ""),
("spindeterminants_psi_det_alpha" , "integer*8" , "(N_int*bit_kind/8,det_alpha_num)"),
("spindeterminants_psi_det_beta" , "integer*8" , "(N_int*bit_kind/8,det_beta_num)"),
("spindeterminants_psi_coef_matrix_rows" , "integer" , "(det_num_input)"),
("spindeterminants_psi_coef_matrix_columns" , "integer" , "(det_num_input)"),
("spindeterminants_psi_coef_matrix_values" , "double precision" , "(det_num_input,N_states)"),
]
def do_subst(t0,d):
t = t0
t = t.replace("$X",d[0])
t = t.replace("$T",d[1])
t = t.replace("$D",d[2])
if d[1].startswith("character"):
size = d[1].split("*")[1][1:-1]
u = "character"
elif d[1].startswith("double precision"):
u = d[1].replace(" ","_")
size = "1"
elif "*" in d[1]:
size = "1"
u = d[1].replace("*","")
else:
size = "1"
u = d[1]
t = t.replace("$U",u)
if d[2] == "":
t = t.replace("$S",size)
else:
if size == "1":
t = t.replace("$S","size(res)")
else:
t = t.replace("$S","%s*size(res)"%(size))
provide = ""
tmp = d[2].replace('(','').replace(')','')
for i in "+-*/":
tmp = tmp.replace(i,',')
for i in tmp.split(','):
if ":" in i:
i = i.split(':')[1]
try:
eval(i+"+1")
except NameError:
provide += " PROVIDE "+i+"\n"
t = t.replace("$P",provide)
print t
t0 = """
subroutine get_$X(res)
implicit none
BEGIN_DOC
! Calls EZFIO subroutine to get $X
END_DOC
$T :: res$D
integer :: ierr, i
logical :: exists
PROVIDE ezfio_filename
$P
if (.not.is_worker) then
call ezfio_has_$X(exists)
if (exists) then
call ezfio_get_$X(res)
call ezfio_free_$X
else
call ezfio_set_$X(res)
endif
else
call zmq_ezfio_has('$X',exists)
if (exists) then
call zmq_ezfio_get_$U('$X',res,$S)
endif
endif
end
"""
t1 = """
subroutine get_$X(res)
implicit none
BEGIN_DOC
! Calls EZFIO subroutine to get $X
END_DOC
$T :: res$D
integer :: ierr
PROVIDE ezfio_filename
$P
if (.not.is_worker) then
call ezfio_get_$X(res)
call ezfio_free_$X
else
call zmq_ezfio_get_$U('$X',res,$S)
endif
end
"""
for i,d in enumerate(data):
do_subst(t0,d)
for i,d in enumerate(data_no_set):
i += len(data)
do_subst(t1,d)
END_SHELL

22
src/finish.irp.f Normal file
View File

@ -0,0 +1,22 @@
subroutine abrt (here,message)
implicit none
character*(*) :: here
character*(*) :: message
write(0,*) '-------------------------'
write(0,*) 'Error in '//trim(here)//':'
write(0,*) '-------------------------'
write(0,*) trim(message)//'.'
write(0,*) '-------------------------'
if (is_worker) then
call worker_log(here,message)
call sleep(2)
endif
call finish()
stop
end
subroutine finish()
implicit none
call ezfio_finish()
end

693
src/mo.irp.f Normal file
View File

@ -0,0 +1,693 @@
BEGIN_PROVIDER [ integer, mo_num ]
&BEGIN_PROVIDER [ integer, mo_num_8 ]
implicit none
BEGIN_DOC
! Number of Molecular orbitals
END_DOC
integer, external :: mod_align
mo_num = maxval(present_mos)
call iinfo(irp_here,'mo_num',mo_num)
mo_num_8 = mod_align(mo_num)
END_PROVIDER
BEGIN_PROVIDER [ real, mo_coef_input, (ao_num_8,mo_tot_num) ]
implicit none
BEGIN_DOC
! Molecular orbital coefficients read from the input file
END_DOC
integer :: i, j
real,allocatable :: buffer(:,:)
allocate (buffer(ao_num,mo_tot_num))
buffer = 0.
call get_mo_basis_mo_coef(buffer)
do i=1,mo_tot_num
do j=1,ao_num
mo_coef_input(j,i) = buffer(j,i)
enddo
do j=ao_num+1,ao_num_8
mo_coef_input(j,i) = 0.
enddo
call set_order(mo_coef_input(1,i),ao_nucl_sort_idx,ao_num)
enddo
deallocate(buffer)
END_PROVIDER
BEGIN_PROVIDER [ real, mo_scale ]
&BEGIN_PROVIDER [ real, mo_norm ]
implicit none
BEGIN_DOC
! Scaling factor for MOs to keep the determinant in a defined domain
END_DOC
mo_scale = 1.d0/(0.4d0*log(float(elec_num+1)))
mo_norm = mo_scale*mo_scale
END_PROVIDER
BEGIN_PROVIDER [ real, mo_coef, (ao_num_8,mo_num_8) ]
implicit none
BEGIN_DOC
! Molecular orbital coefficients
END_DOC
integer :: i, j
do j=1,mo_num
do i=1,ao_num_8
mo_coef(i,j) = mo_coef_input(i,j)
enddo
enddo
do j =mo_num+1,mo_num_8
!DIR$ VECTOR ALIGNED
do i=1,ao_num_8
mo_coef(i,j) = 0.
enddo
enddo
! Input MOs are not needed any more
FREE mo_coef_input
real :: f
f = 1./mo_scale
do j=1,mo_num
!DIR$ VECTOR ALIGNED
!DIR$ LOOP COUNT (2000)
do i=1,ao_num_8
mo_coef(i,j) *= f
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ real, mo_coef_transp, (mo_num_8,ao_num_8) ]
implicit none
BEGIN_DOC
! Transpose of the Molecular orbital coefficients
END_DOC
call transpose(mo_coef,ao_num_8,mo_coef_transp,mo_num_8,ao_num_8,mo_num_8)
END_PROVIDER
BEGIN_PROVIDER [ integer, mo_coef_transp_non_zero_idx, (0:mo_num,ao_num) ]
&BEGIN_PROVIDER [ real, mo_coef_transp_sparsity ]
implicit none
BEGIN_DOC
! Indices of the non-zero elements of the transpose of the Molecular
! orbital coefficients
END_DOC
integer :: i, j
integer :: idx
mo_coef_transp_sparsity = 0.
do j=1,ao_num
idx = 0
do i=1,mo_num
if (mo_coef_transp(i,j) /= 0.) then
idx += 1
mo_coef_transp_non_zero_idx(idx,j) = i
endif
enddo
mo_coef_transp_non_zero_idx(0,j) = idx
mo_coef_transp_sparsity += float(idx)
enddo
mo_coef_transp_sparsity *= 1./(mo_num*ao_num)
END_PROVIDER
BEGIN_PROVIDER [ real, mo_coef_transp_present, (num_present_mos_8,ao_num_8) ]
implicit none
BEGIN_DOC
! mo_coef_transp without MOs absent in all determinants
END_DOC
integer :: i,j,n
mo_coef_transp_present = 0.
do i=1,ao_num
do j=1,num_present_mos
mo_coef_transp_present(j,i) = mo_coef_transp(present_mos(j),i)
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ real, mo_value_transp, ((-simd_sp+1):mo_num_8,elec_num) ]
&BEGIN_PROVIDER [ real, mo_grad_transp_x, ((-2*simd_sp+1):mo_num_8,elec_num) ]
&BEGIN_PROVIDER [ real, mo_grad_transp_y, ((-3*simd_sp+1):mo_num_8,elec_num) ]
&BEGIN_PROVIDER [ real, mo_grad_transp_z, ((-4*simd_sp+1):mo_num_8,elec_num) ]
&BEGIN_PROVIDER [ real, mo_lapl_transp, ((-5*simd_sp+1):mo_num_8,elec_num) ]
implicit none
BEGIN_DOC
! Values, gradients, laplacians of the molecular orbitals
!
! Arrays are padded for efficiency
END_DOC
integer :: i, j, k, l, m
PROVIDE primitives_reduced
if (do_nucl_fitcusp) then
PROVIDE nucl_fitcusp_param
PROVIDE nucl_elec_dist_vec
PROVIDE nucl_elec_dist_inv
endif
do i=1,elec_num
if (i>1) then
ao_elec = i
TOUCH ao_elec
endif
if (num_present_mos == mo_num) then
call sparse_full_mv(mo_coef_transp,mo_num_8, &
ao_value_block(1),ao_num_8, &
ao_grad_block_x(1), &
ao_grad_block_y(1), &
ao_grad_block_z(1), &
ao_lapl_block(1), &
ao_value_non_zero_idx(0), &
mo_value_transp(1,i),mo_num_8, &
mo_grad_transp_x(1,i), &
mo_grad_transp_y(1,i), &
mo_grad_transp_z(1,i), &
mo_lapl_transp(1,i), &
ao_num)
else
call sparse_full_mv(mo_coef_transp_present,num_present_mos_8, &
ao_value_block(1),ao_num_8, &
ao_grad_block_x(1), &
ao_grad_block_y(1), &
ao_grad_block_z(1), &
ao_lapl_block(1), &
ao_value_non_zero_idx(0), &
mo_value_transp(1,i),mo_num_8, &
mo_grad_transp_x(1,i), &
mo_grad_transp_y(1,i), &
mo_grad_transp_z(1,i), &
mo_lapl_transp(1,i), &
ao_num)
do j=num_present_mos,1,-1
mo_value_transp (present_mos(j),i) = mo_value_transp (j,i)
mo_grad_transp_x(present_mos(j),i) = mo_grad_transp_x(j,i)
mo_grad_transp_y(present_mos(j),i) = mo_grad_transp_y(j,i)
mo_grad_transp_z(present_mos(j),i) = mo_grad_transp_z(j,i)
mo_lapl_transp (present_mos(j),i) = mo_lapl_transp (j,i)
if (present_mos(j) == j) then
exit
endif
enddo
endif
if (do_nucl_fitcusp) then
real :: r, r2, r_inv, d, expzr, Z, Z2, a, b, c, phi, rx, ry, rz
do k=1,nucl_num
r = nucl_elec_dist(k,i)
if (r > nucl_fitcusp_radius(k)) then
cycle
endif
r_inv = nucl_elec_dist_inv(k,i)
!DIR$ LOOP COUNT (500)
do j=1,mo_num
mo_value_transp(j,i) = mo_value_transp(j,i) + nucl_fitcusp_param(1,j,k) +&
r * (nucl_fitcusp_param(2,j,k) + &
r * (nucl_fitcusp_param(3,j,k) + &
r * nucl_fitcusp_param(4,j,k) ))
mo_lapl_transp(j,i) = mo_lapl_transp(j,i) + &
nucl_fitcusp_param(2,j,k)*(r_inv+r_inv) + &
6.*nucl_fitcusp_param(3,j,k) + &
r * 12.*nucl_fitcusp_param(4,j,k)
c = r_inv * (nucl_fitcusp_param(2,j,k) + &
r * (2.*nucl_fitcusp_param(3,j,k) + &
r * 3.*nucl_fitcusp_param(4,j,k) ))
mo_grad_transp_x(j,i) = mo_grad_transp_x(j,i) + nucl_elec_dist_vec(1,k,i)*c
mo_grad_transp_y(j,i) = mo_grad_transp_y(j,i) + nucl_elec_dist_vec(2,k,i)*c
mo_grad_transp_z(j,i) = mo_grad_transp_z(j,i) + nucl_elec_dist_vec(3,k,i)*c
enddo
exit
enddo ! k
endif
enddo ! i
ao_elec = 1
SOFT_TOUCH ao_elec
if (do_prepare) then
real :: lambda, t
! Scale off-diagonal elements
t = prepare_walkers_t
do i=1,mo_num
!DIR$ LOOP COUNT (100)
do j=1,elec_alpha_num
if (i /= j) then
mo_value_transp(i,j) *= t
mo_grad_transp_x(i,j) *= t
mo_grad_transp_y(i,j) *= t
mo_grad_transp_z(i,j) *= t
mo_lapl_transp(i,j) *= t
endif
enddo
!DIR$ LOOP COUNT (100)
do j=1,elec_beta_num
if (i /= j) then
mo_value_transp(i,j+elec_alpha_num) *= t
mo_grad_transp_x(i,j+elec_alpha_num) *= t
mo_grad_transp_y(i,j+elec_alpha_num) *= t
mo_grad_transp_z(i,j+elec_alpha_num) *= t
mo_lapl_transp(i,j+elec_alpha_num) *= t
endif
enddo
enddo
endif
END_PROVIDER
BEGIN_PROVIDER [ real, mo_value, (elec_num_8,mo_num) ]
implicit none
BEGIN_DOC
! Values of the molecular orbitals
END_DOC
integer :: i,j
integer, save :: ifirst = 0
if (ifirst == 0) then
ifirst = 1
PROVIDE primitives_reduced
!DIR$ VECTOR ALIGNED
mo_value = 0.
endif
call transpose(mo_value_transp(1,1),mo_num_8+simd_sp,mo_value,elec_num_8,mo_num,elec_num)
END_PROVIDER
BEGIN_PROVIDER [ double precision, mo_grad_x, (elec_num_8,mo_num) ]
&BEGIN_PROVIDER [ double precision, mo_grad_y, (elec_num_8,mo_num) ]
&BEGIN_PROVIDER [ double precision, mo_grad_z, (elec_num_8,mo_num) ]
implicit none
BEGIN_DOC
! Gradients of the molecular orbitals
END_DOC
integer :: i,j
integer, save :: ifirst = 0
if (ifirst == 0) then
!DIR$ VECTOR ALIGNED
mo_grad_x = 0.d0
!DIR$ VECTOR ALIGNED
mo_grad_y = 0.d0
!DIR$ VECTOR ALIGNED
mo_grad_z = 0.d0
ifirst = 1
PROVIDE primitives_reduced
endif
! Transpose x last for cache efficiency
call transpose_to_dp(mo_grad_transp_y(1,1),mo_num_8+3*simd_sp,mo_grad_y(1,1),elec_num_8,mo_num,elec_num)
call transpose_to_dp(mo_grad_transp_z(1,1),mo_num_8+4*simd_sp,mo_grad_z(1,1),elec_num_8,mo_num,elec_num)
call transpose_to_dp(mo_grad_transp_x(1,1),mo_num_8+2*simd_sp,mo_grad_x(1,1),elec_num_8,mo_num,elec_num)
END_PROVIDER
BEGIN_PROVIDER [ double precision, mo_grad_lapl, (4,elec_num,mo_num) ]
implicit none
BEGIN_DOC
! Gradients and laplacian
END_DOC
integer :: i,j
do j=1,mo_num
do i=1,elec_num
mo_grad_lapl(1,i,j) = mo_grad_x(i,j)
mo_grad_lapl(2,i,j) = mo_grad_y(i,j)
mo_grad_lapl(3,i,j) = mo_grad_z(i,j)
mo_grad_lapl(4,i,j) = mo_lapl (i,j)
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, mo_lapl, (elec_num_8,mo_num) ]
implicit none
BEGIN_DOC
! Laplacians of the molecular orbitals
END_DOC
integer :: i,j
integer, save :: ifirst = 0
if (ifirst == 0) then
ifirst = 1
PROVIDE primitives_reduced
!DIR$ VECTOR ALIGNED
mo_lapl = 0.d0
endif
call transpose_to_dp(mo_lapl_transp(1,1),mo_num_8+5*simd_sp,mo_lapl,elec_num_8,mo_num,elec_num)
END_PROVIDER
BEGIN_PROVIDER [ real, prepare_walkers_t ]
implicit none
BEGIN_DOC
! prepare_walkers_t : scaling of the off-diagonal elements
! of the mo_value matrix
END_DOC
prepare_walkers_t = 1.
END_PROVIDER
BEGIN_PROVIDER [ integer, mo_tot_num ]
BEGIN_DOC
! Total number of MOs in the EZFIO file
END_DOC
mo_tot_num = -1
call get_mo_basis_mo_tot_num(mo_tot_num)
if (mo_tot_num <= 0) then
call abrt(irp_here,'Total number of MOs can''t be <0')
endif
call iinfo(irp_here,'mo_tot_num',mo_tot_num)
END_PROVIDER
!-----------------
! Fit cusp
!-----------------
BEGIN_PROVIDER [ double precision , mo_value_at_nucl, (mo_num_8,nucl_num) ]
implicit none
BEGIN_DOC
! Values of the molecular orbitals at the nucleus without the
! S components of the current nucleus
END_DOC
integer :: i, j, k, l
real :: ao_value_at_nucl_no_S(ao_num)
do k=1,nucl_num
point(1) = nucl_coord(k,1)
point(2) = nucl_coord(k,2)
point(3) = nucl_coord(k,3)
TOUCH point
PROVIDE ao_value_p
!DIR$ LOOP COUNT (2000)
do i=1,ao_num
if (ao_nucl(i) /= k) then
ao_value_at_nucl_no_S(i) = ao_value_p(i)
else
ao_value_at_nucl_no_S(i) = 0.
endif
enddo
integer :: jj
do jj=1,num_present_mos
j = present_mos(jj)
mo_value_at_nucl(j,k) = 0.
!DIR$ VECTOR ALIGNED
do i=1,ao_num
mo_value_at_nucl(j,k) = mo_value_at_nucl(j,k) + mo_coef(i,j)*ao_value_at_nucl_no_S(i)
enddo
enddo
enddo
FREE ao_value_p ao_grad_p ao_lapl_p ao_axis_grad_p ao_oned_grad_p ao_oned_prim_grad_p ao_oned_lapl_p ao_axis_lapl_p ao_oned_prim_lapl_p ao_oned_p ao_oned_prim_p ao_axis_p ao_axis_power_p
SOFT_TOUCH point
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_value_at_fitcusp_radius, (ao_num_8,nucl_num) ]
&BEGIN_PROVIDER [ double precision, ao_grad_at_fitcusp_radius, (ao_num_8,nucl_num) ]
&BEGIN_PROVIDER [ double precision, ao_lapl_at_fitcusp_radius, (ao_num_8,nucl_num) ]
implicit none
BEGIN_DOC
! Values of the atomic orbitals without S components on atoms
END_DOC
integer :: i, j, k
do k=1,nucl_num
point(1) = nucl_coord(k,1)
point(2) = nucl_coord(k,2)
point(3) = nucl_coord(k,3)+ nucl_fitcusp_radius(k)
TOUCH point
!DIR$ LOOP COUNT (2000)
do j=1,ao_num
ao_value_at_fitcusp_radius(j,k) = ao_value_p(j)
ao_grad_at_fitcusp_radius(j,k) = ao_grad_p(j,3)
ao_lapl_at_fitcusp_radius(j,k) = ao_lapl_p(j)
if ( (ao_nucl(j) /= k).or.(ao_power(j,4) >0) ) then
ao_value_at_fitcusp_radius(j,k) = 0.
ao_grad_at_fitcusp_radius(j,k) = 0.
ao_lapl_at_fitcusp_radius(j,k) = 0.
endif
enddo
enddo
FREE ao_value_p ao_grad_p ao_lapl_p ao_axis_grad_p ao_oned_grad_p ao_oned_prim_grad_p ao_oned_lapl_p ao_axis_lapl_p ao_oned_prim_lapl_p ao_oned_p ao_oned_prim_p ao_axis_p ao_axis_power_p
SOFT_TOUCH point
END_PROVIDER
BEGIN_PROVIDER [ double precision, mo_value_at_fitcusp_radius, (mo_num_8,nucl_num) ]
&BEGIN_PROVIDER [ double precision, mo_grad_at_fitcusp_radius, (mo_num_8,nucl_num) ]
&BEGIN_PROVIDER [ double precision, mo_lapl_at_fitcusp_radius, (mo_num_8,nucl_num) ]
implicit none
BEGIN_DOC
! Values of the molecular orbitals without S components on atoms
END_DOC
integer :: i, j, k, l
do k=1,nucl_num
do j=1,mo_num
mo_value_at_fitcusp_radius(j,k) = 0.d0
mo_grad_at_fitcusp_radius(j,k) = 0.d0
mo_lapl_at_fitcusp_radius(j,k) = 0.d0
!DIR$ VECTOR ALIGNED
do i=1,ao_num
mo_value_at_fitcusp_radius(j,k) = mo_value_at_fitcusp_radius(j,k) + mo_coef(i,j)*ao_value_at_fitcusp_radius(i,k)
mo_grad_at_fitcusp_radius(j,k) = mo_grad_at_fitcusp_radius(j,k) + mo_coef(i,j)*ao_grad_at_fitcusp_radius(i,k)
mo_lapl_at_fitcusp_radius(j,k) = mo_lapl_at_fitcusp_radius(j,k) + mo_coef(i,j)*ao_lapl_at_fitcusp_radius(i,k)
enddo
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ real, nucl_fitcusp_param, (4,mo_num,nucl_num) ]
implicit none
BEGIN_DOC
! Parameters of the splines
END_DOC
integer :: i,k, niter
character*(80) :: message
nucl_fitcusp_param = 0.d0
do k=1,nucl_num
double precision :: r, Z
Z = nucl_charge(k)
if (Z < 1.d-2) then
! Avoid dummy atoms
cycle
endif
R = nucl_fitcusp_radius(k)
!DIR$ LOOP COUNT (500)
do i=1,mo_num
double precision :: lap_phi, grad_phi, phi, eta
lap_phi = mo_lapl_at_fitcusp_radius(i,k)
grad_phi = mo_grad_at_fitcusp_radius(i,k)
phi = mo_value_at_fitcusp_radius(i,k)
eta = mo_value_at_nucl(i,k)
nucl_fitcusp_param(1,i,k) = -(R*(2.d0*eta*Z-6.d0*grad_phi)+lap_phi*R*R+6.d0*phi)/(2.d0*R*Z-6.d0)
nucl_fitcusp_param(2,i,k) = (lap_phi*R*R*Z-6.d0*grad_phi*R*Z+6.d0*phi*Z+6.d0*eta*Z)/(2.d0*R*Z-6.d0)
nucl_fitcusp_param(3,i,k) = -(R*(-5.d0*grad_phi*Z-1.5d0*lap_phi)+lap_phi*R*R*Z+3.d0*phi*Z+&
3.d0*eta*Z+6.d0*grad_phi)/(R*R*Z-3.d0*R)
nucl_fitcusp_param(4,i,k) = (R*(-2.d0*grad_phi*Z-lap_phi)+0.5d0*lap_phi*R*R*Z+phi*Z+&
eta*Z+3.d0*grad_phi)/(R*R*R*Z-3.d0*R*R)
enddo
enddo
END_PROVIDER
subroutine sparse_full_mv(A,LDA, &
B1,LDB, &
B2, B3, B4, B5, indices, &
C1,LDC,C2,C3,C4,C5,an)
implicit none
BEGIN_DOC
! Performs a vectorized product between a dense matrix (the MO coefficients
! matrix) and 5 sparse vectors (the value, gradients and laplacian of the AOs).
END_DOC
integer, intent(in) :: an,LDA,LDB,LDC
integer, intent(in) :: indices(0:LDB)
real, intent(in) :: A(LDA,an)
real, intent(in) :: B1(LDB)
real, intent(in) :: B2(LDB)
real, intent(in) :: B3(LDB)
real, intent(in) :: B4(LDB)
real, intent(in) :: B5(LDB)
real, intent(out) :: C1(LDC)
real, intent(out) :: C2(LDC)
real, intent(out) :: C3(LDC)
real, intent(out) :: C4(LDC)
real, intent(out) :: C5(LDC)
!DIR$ ASSUME_ALIGNED A : $IRP_ALIGN
!DIR$ ASSUME_ALIGNED B1 : $IRP_ALIGN
!DIR$ ASSUME_ALIGNED B2 : $IRP_ALIGN
!DIR$ ASSUME_ALIGNED B3 : $IRP_ALIGN
!DIR$ ASSUME_ALIGNED B4 : $IRP_ALIGN
!DIR$ ASSUME_ALIGNED B5 : $IRP_ALIGN
!DIR$ ASSUME_ALIGNED C1 : $IRP_ALIGN
!DIR$ ASSUME_ALIGNED C2 : $IRP_ALIGN
!DIR$ ASSUME_ALIGNED C3 : $IRP_ALIGN
!DIR$ ASSUME_ALIGNED C4 : $IRP_ALIGN
!DIR$ ASSUME_ALIGNED C5 : $IRP_ALIGN
integer :: kao, kmax, kmax2, kmax3
integer :: i,j,k
integer :: k_vec(8)
!DIR$ ATTRIBUTES ALIGN: $IRP_ALIGN :: k_vec
real :: d11, d12, d13, d14, d15
real :: d21, d22, d23, d24, d25
real :: d31, d32, d33, d34, d35
real :: d41, d42, d43, d44, d45
! LDC and LDA have to be factors of simd_sp
!DIR$ VECTOR ALIGNED
!DIR$ LOOP COUNT (256)
!$OMP SIMD
do j=1,LDC
C1(j) = 0.
C2(j) = 0.
C3(j) = 0.
C4(j) = 0.
C5(j) = 0.
enddo
!$OMP END SIMD
kmax2 = ishft(indices(0),-2)
kmax2 = ishft(kmax2,2)
kmax3 = indices(0)
!DIR$ LOOP COUNT (200)
do kao=1,kmax2,4
k_vec(1) = indices(kao )
k_vec(2) = indices(kao+1)
k_vec(3) = indices(kao+2)
k_vec(4) = indices(kao+3)
d11 = B1(kao )
d21 = B1(kao+1)
d31 = B1(kao+2)
d41 = B1(kao+3)
d12 = B2(kao )
d22 = B2(kao+1)
d32 = B2(kao+2)
d42 = B2(kao+3)
!DIR$ LOOP COUNT (256)
do k=0,LDA-1,$IRP_ALIGN/4
!DIR$ VECTOR ALIGNED
!$OMP SIMD
do j=1,$IRP_ALIGN/4
C1(j+k) = C1(j+k) + A(j+k,k_vec(1))*d11 + A(j+k,k_vec(2))*d21&
+ A(j+k,k_vec(3))*d31 + A(j+k,k_vec(4))*d41
C2(j+k) = C2(j+k) + A(j+k,k_vec(1))*d12 + A(j+k,k_vec(2))*d22&
+ A(j+k,k_vec(3))*d32 + A(j+k,k_vec(4))*d42
enddo
!$OMP END SIMD
enddo
d13 = B3(kao )
d23 = B3(kao+1)
d33 = B3(kao+2)
d43 = B3(kao+3)
d14 = B4(kao )
d24 = B4(kao+1)
d34 = B4(kao+2)
d44 = B4(kao+3)
!DIR$ LOOP COUNT (256)
do k=0,LDA-1,$IRP_ALIGN/4
!DIR$ VECTOR ALIGNED
!$OMP SIMD
do j=1,$IRP_ALIGN/4
C3(j+k) = C3(j+k) + A(j+k,k_vec(1))*d13 + A(j+k,k_vec(2))*d23&
+ A(j+k,k_vec(3))*d33 + A(j+k,k_vec(4))*d43
C4(j+k) = C4(j+k) + A(j+k,k_vec(1))*d14 + A(j+k,k_vec(2))*d24&
+ A(j+k,k_vec(3))*d34 + A(j+k,k_vec(4))*d44
enddo
!$OMP END SIMD
enddo
d15 = B5(kao )
d25 = B5(kao+1)
d35 = B5(kao+2)
d45 = B5(kao+3)
!DIR$ LOOP COUNT (256)
do k=0,LDA-1,$IRP_ALIGN/4
!DIR$ VECTOR ALIGNED
!$OMP SIMD
do j=1,$IRP_ALIGN/4
C5(j+k) = C5(j+k) + A(j+k,k_vec(1))*d15 + A(j+k,k_vec(2))*d25&
+ A(j+k,k_vec(3))*d35 + A(j+k,k_vec(4))*d45
enddo
!$OMP END SIMD
enddo
enddo
!DIR$ LOOP COUNT (200)
do kao = kmax2+1, kmax3
k_vec(1) = indices(kao)
d11 = B1(kao)
d12 = B2(kao)
d13 = B3(kao)
d14 = B4(kao)
d15 = B5(kao)
!DIR$ VECTOR ALIGNED
!DIR$ LOOP COUNT (256)
do k=0,LDA-1,simd_sp
!DIR$ VECTOR ALIGNED
!$OMP SIMD
do j=1,$IRP_ALIGN/4
C1(j+k) = C1(j+k) + A(j+k,k_vec(1))*d11
C2(j+k) = C2(j+k) + A(j+k,k_vec(1))*d12
C3(j+k) = C3(j+k) + A(j+k,k_vec(1))*d13
C4(j+k) = C4(j+k) + A(j+k,k_vec(1))*d14
C5(j+k) = C5(j+k) + A(j+k,k_vec(1))*d15
enddo
!$OMP END SIMD
enddo
enddo
end

21
src/mo_point.irp.f Normal file
View File

@ -0,0 +1,21 @@
BEGIN_PROVIDER [ real, mo_value_p, (mo_tot_num) ]
implicit none
BEGIN_DOC
! Values of the molecular orbitals
END_DOC
integer :: i, j, k
do j=1,mo_num
mo_value_p(j) = 0.
enddo
do k=1,ao_num
do j=1,mo_num
mo_value_p(j) = mo_value_p(j)+mo_coef_transp(j,k)*ao_value_p(k)
enddo
enddo
END_PROVIDER

148
src/nuclei.irp.f Normal file
View File

@ -0,0 +1,148 @@
BEGIN_PROVIDER [ integer, nucl_num ]
&BEGIN_PROVIDER [ integer, nucl_num_8 ]
implicit none
BEGIN_DOC
! Number of nuclei
END_DOC
nucl_num = -1
call get_nuclei_nucl_num(nucl_num)
if (nucl_num <= 0) then
call abrt(irp_here,'Number of nuclei should be > 0')
endif
integer, external :: mod_align
nucl_num_8 = mod_align(nucl_num)
END_PROVIDER
BEGIN_PROVIDER [ real, nucl_charge, (nucl_num) ]
implicit none
BEGIN_DOC
! Nuclear charge
END_DOC
nucl_charge = -1.d0
call get_nuclei_nucl_charge(nucl_charge)
integer :: i
do i=1,nucl_num
if (nucl_charge(i) < 0.) then
call abrt(irp_here,'Nuclear charges should be > 0')
endif
enddo
END_PROVIDER
BEGIN_PROVIDER [ real, nucl_coord, (nucl_num_8,3) ]
implicit none
BEGIN_DOC
! Nuclear coordinates
END_DOC
nucl_coord = 0.
real, allocatable :: buffer(:,:)
allocate (buffer(nucl_num,3))
buffer = 0.
call get_nuclei_nucl_coord(buffer)
integer :: i,j
do i=1,3
do j=1,nucl_num
nucl_coord(j,i) = buffer(j,i)
enddo
enddo
deallocate(buffer)
END_PROVIDER
BEGIN_PROVIDER [ real, nucl_coord_transp, (8,nucl_num)
implicit none
BEGIN_DOC
! Transposed array of nucl_coord
END_DOC
integer :: i, k
integer, save :: ifirst = 0
if (ifirst == 0) then
ifirst = 1
nucl_coord_transp = 0.
endif
!DIR$ VECTOR ALIGNED
do i=1,nucl_num
nucl_coord_transp(1,i) = nucl_coord(i,1)
nucl_coord_transp(2,i) = nucl_coord(i,2)
nucl_coord_transp(3,i) = nucl_coord(i,3)
enddo
END_PROVIDER
BEGIN_PROVIDER [ real, nucl_dist, (nucl_num_8,nucl_num) ]
&BEGIN_PROVIDER [ real, nucl_dist_vec_x, (nucl_num_8,nucl_num) ]
&BEGIN_PROVIDER [ real, nucl_dist_vec_y, (nucl_num_8,nucl_num) ]
&BEGIN_PROVIDER [ real, nucl_dist_vec_z, (nucl_num_8,nucl_num) ]
implicit none
BEGIN_DOC
! nucl_dist : Nucleus-nucleus distances : nucl_dist(i,j) = |R_i-R_j|
! nucl_dist_vec : Nucleus-nucleus distances vectors
END_DOC
integer :: ie1, ie2, l
integer,save :: ifirst = 0
if (ifirst == 0) then
ifirst = 1
nucl_dist = 0.
nucl_dist_vec_x = 0.
nucl_dist_vec_y = 0.
nucl_dist_vec_z = 0.
endif
do ie2 = 1,nucl_num
!DIR$ VECTOR ALIGNED
!DIR$ LOOP COUNT (100)
do ie1 = 1,nucl_num
nucl_dist_vec_x(ie1,ie2) = nucl_coord(ie1,1) - nucl_coord(ie2,1)
nucl_dist_vec_y(ie1,ie2) = nucl_coord(ie1,2) - nucl_coord(ie2,2)
nucl_dist_vec_z(ie1,ie2) = nucl_coord(ie1,3) - nucl_coord(ie2,3)
enddo
!DIR$ VECTOR ALIGNED
!DIR$ LOOP COUNT (100)
do ie1 = 1,nucl_num
nucl_dist (ie1,ie2) = nucl_dist_vec_x(ie1,ie2)*nucl_dist_vec_x(ie1,ie2) +&
nucl_dist_vec_y(ie1,ie2)*nucl_dist_vec_y(ie1,ie2) + &
nucl_dist_vec_z(ie1,ie2)*nucl_dist_vec_z(ie1,ie2)
nucl_dist(ie1,ie2) = sqrt(nucl_dist (ie1,ie2))
ASSERT (nucl_dist(ie1,ie2) > 0.)
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ real, nucl_fitcusp_radius, (nucl_num) ]
implicit none
BEGIN_DOC
! Distance threshold for the fit
END_DOC
real :: def(nucl_num)
integer :: k
if (.not.do_nucl_fitcusp) then
nucl_fitcusp_radius = 0.d0
return
endif
do k=1,nucl_num
nucl_fitcusp_radius(k) = .5/nucl_charge(k)
enddo
call get_nuclei_nucl_fitcusp_radius(nucl_fitcusp_radius)
! Avoid dummy atoms
do k=1,nucl_num
if (nucl_charge(k) < 5.d-1) then
nucl_fitcusp_radius(k) = 0.
endif
enddo
END_PROVIDER

39
src/point.irp.f Normal file
View File

@ -0,0 +1,39 @@
BEGIN_PROVIDER [ real, point, (3) ]
implicit none
BEGIN_DOC
! Coordinates of the current point
END_DOC
point(1) = 0.
point(2) = 0.
point(3) = 0.
END_PROVIDER
BEGIN_PROVIDER [ real, point_nucl_dist_vec, (nucl_num,3) ]
implicit none
BEGIN_DOC
! Distance vector between the current point and the nuclei
END_DOC
integer :: k
do k=1,nucl_num
point_nucl_dist_vec(k,1) = point(1)-nucl_coord(k,1)
point_nucl_dist_vec(k,2) = point(2)-nucl_coord(k,2)
point_nucl_dist_vec(k,3) = point(3)-nucl_coord(k,3)
enddo
END_PROVIDER
BEGIN_PROVIDER [ real, point_nucl_dist, (nucl_num) ]
implicit none
BEGIN_DOC
! Distance between the current point and the nuclei
END_DOC
integer :: k,l
do k=1,nucl_num
point_nucl_dist(k) = point_nucl_dist_vec(k,1)*point_nucl_dist_vec(k,1) +&
point_nucl_dist_vec(k,2)*point_nucl_dist_vec(k,2) + &
point_nucl_dist_vec(k,3)*point_nucl_dist_vec(k,3)
point_nucl_dist(k) = sqrt(point_nucl_dist(k))
enddo
END_PROVIDER

217
src/prepare_walkers.irp.f Normal file
View File

@ -0,0 +1,217 @@
subroutine draw_init_points
implicit none
BEGIN_DOC
! Place randomly electron around nuclei
END_DOC
integer :: iwalk
logical, allocatable :: do_elec(:)
integer :: acc_num
real, allocatable :: xmin(:,:)
integer :: i, j, k, l, kk
real :: norm
allocate (do_elec(elec_num), xmin(3,elec_num))
xmin = -huge(1.)
norm = 0.
do i=1,elec_alpha_num
do j=1,ao_num
norm += mo_coef_transp(i,j)*mo_coef_transp(i,j)
enddo
enddo
norm = sqrt(norm/float(elec_alpha_num))
call rinfo( irp_here, 'Norm : ', norm )
call rinfo( irp_here, 'mo_scale: ' , mo_scale )
mo_coef_transp = mo_coef_transp/norm
double precision :: qmc_ranf
real :: mo_max
do i=1,elec_alpha_num
l=1
xmin(1,i) = mo_coef_transp(i,1)*mo_coef_transp(i,1) - 0.001*qmc_ranf()
do j=2,ao_num
xmin(2,i) = mo_coef_transp(i,j)*mo_coef_transp(i,j) - 0.001*qmc_ranf()
if (xmin(2,i) > xmin(1,i) ) then
xmin(1,i) = xmin(2,i)
l = ao_nucl(j)
endif
enddo
xmin(1,i) = nucl_coord(l,1)
xmin(2,i) = nucl_coord(l,2)
xmin(3,i) = nucl_coord(l,3)
enddo
call iinfo(irp_here, 'Det num = ', det_num )
do k=1,elec_beta_num
i = k+elec_alpha_num
l=1
xmin(1,i) = mo_coef_transp(k,1)*mo_coef_transp(k,1) - 0.001*qmc_ranf()
do j=2,ao_num
xmin(2,i) = mo_coef_transp(k,j)*mo_coef_transp(k,j) - 0.001*qmc_ranf()
if (xmin(2,i) > xmin(1,i) ) then
xmin(1,i) = xmin(2,i)
l = ao_nucl(j)
endif
enddo
xmin(1,i) = nucl_coord(l,1)
xmin(2,i) = nucl_coord(l,2)
xmin(3,i) = nucl_coord(l,3)
enddo
call rinfo( irp_here, 'time step =', time_step )
do iwalk=1,walk_num
call iinfo( irp_here, 'Generate initial positions for walker', iwalk )
acc_num = 0
do_elec = .True.
do while (acc_num < elec_num)
double precision :: gauss
real :: re_compute
re_compute = 0.
do while (re_compute < 1.e-6)
do i=1,elec_num
if (do_elec(i)) then
do l=1,3
elec_coord(i,l) = xmin(l,i) + 2.*(0.5-qmc_ranf())
enddo
endif
enddo
TOUCH elec_coord
re_compute = minval(nucl_elec_dist(1:nucl_num,1:elec_num))
enddo
do i=1,elec_alpha_num
if (do_elec(i)) then
if ( mo_value_transp(i,i)**2 >= qmc_ranf()) then
acc_num += 1
do_elec(i) = .False.
endif
endif
enddo
do i=1,elec_beta_num
if (do_elec(i+elec_alpha_num)) then
if ( mo_value_transp(i,i+elec_alpha_num)**2 >= qmc_ranf()) then
acc_num += 1
do_elec(i+elec_alpha_num) = .False.
endif
endif
enddo
enddo
do l=1,3
do i=1,elec_num+1
elec_coord_full(i,l,iwalk) = elec_coord(i,l)
enddo
enddo
enddo
if (.not.is_worker) then
call ezfio_set_electrons_elec_coord_pool_size(walk_num)
call ezfio_set_electrons_elec_coord_pool(elec_coord_full)
endif
SOFT_TOUCH elec_coord elec_coord_full
deallocate (do_elec, xmin)
end
subroutine run_prepare_walkers
implicit none
BEGIN_DOC
! Create starting points for walkers
END_DOC
include 'types.F'
integer :: istep, iwalk
integer :: i,j, l
do iwalk=1,walk_num
do l=1,3
do i=1,elec_num+1
elec_coord(i,l) = elec_coord_full(i,l,iwalk)
enddo
enddo
TOUCH elec_coord
double precision :: qmc_ranf, rcond, lambda
rcond = 100.d0
lambda = 1.d0
do while ( (rcond > 3.d0) .or. (rcond < -3.d0) )
rcond = 0.
do i=1,elec_alpha_num
rcond += log(lambda*abs(mo_value_transp(i,i)))
enddo
do i=1,elec_beta_num
rcond += log(lambda*abs(mo_value_transp(i,elec_alpha_num+i)))
enddo
if (rcond > 2.d0) then
lambda = lambda/(1.d0+.1*qmc_ranf())
endif
if (rcond< -2.d0) then
lambda = lambda*(1.d0+.1*qmc_ranf())
endif
enddo
do i=1,ao_num
!DIR$ VECTOR ALIGNED
do j=1,mo_num_8
mo_coef_transp(j,i) *= lambda
enddo
enddo
TOUCH mo_coef_transp
call iinfo (irp_here, 'Starting walker ', iwalk )
do istep=1,1000
if (single_det_value == 0.d0) then
exit
endif
prepare_walkers_t = float(istep)/1000.
TOUCH prepare_walkers_t
rcond = log(abs(dble(single_det_value)))
real :: factor
rcond = log(abs(dble(single_det_value)))
integer :: icount
icount = 0
do while ( (rcond > 10.d0) .or. (rcond < -10.d0) )
icount += 1
if (icount == 1000) then
exit
endif
if (rcond > 10.d0) then
factor = 1./(1.+.10) !*qmc_ranf())
else if (rcond< -10.d0) then
factor = 1.+.10 !*qmc_ranf()
endif
do j=1,ao_num
!DIR$ VECTOR ALIGNED
do i=1,mo_num_8
mo_coef_transp(i,j) *= factor
enddo
enddo
TOUCH mo_coef_transp
rcond = log(abs(dble(single_det_value)))
enddo
double precision :: p,q
logical :: accepted
real :: delta_x
accepted = .False.
do while (.not.accepted)
if (vmc_algo == t_Brownian) then
call brownian_step(p,q,accepted,delta_x)
else if (vmc_algo == t_Langevin) then
call langevin_step(p,q,accepted,delta_x)
endif
enddo
enddo
do l=1,3
do i=1,elec_num+1
elec_coord_full(i,l,iwalk) = elec_coord(i,l)
enddo
enddo
call iinfo(irp_here, 'Walker done', iwalk)
TOUCH elec_coord_full
enddo
end

95
src/properties.py Executable file
View File

@ -0,0 +1,95 @@
#!/usr/bin/python
import string
import os
properties = []
dims = {}
files = filter(lambda x: x.startswith("PROPERTIES") and \
x.endswith("irp.f"), os.listdir(os.getcwd()))
files = map(lambda x: 'PROPERTIES/'+x, filter(lambda x: x.endswith("irp.f"), os.listdir(os.getcwd()+'/PROPERTIES')))
#files = filter(lambda x: x.endswith("irp.f"), os.listdir(os.getcwd()))
for filename in files:
lines = []
check_dims = False
file = open(filename,'r')
lines += file.readlines()
file.close()
for i,line in enumerate(lines):
if line.startswith("! PROPERTIES"):
lines = lines[i:]
break
for line in map(lambda x: x.lower(),lines):
if line.lstrip().startswith('begin_provider'):
check_dims = False
buffer = line
buffer = buffer.split('[')[1]
buffer = buffer.split(']')[0]
buffer = buffer.split(',')
if (len(buffer) == 2):
buffer.append("")
else:
buffer = [ buffer[0], buffer[1], ','.join(buffer[2:]) ]
check_dims = True
buffer = map(string.strip,buffer)
properties.append(buffer)
current_prop = buffer[1]
elif check_dims:
if 'dimensions :' in line:
dims[current_prop] = line.split(':')[1].strip()
def sq(item):
return [item[0], item[1]+"_2", item[2]]
properties_with_square = properties + map(sq,properties)
for p in [ properties, properties_with_square ]:
def c(x,y):
if x[1] > y[1]: return 1
if x[1] == y[1]: return 0
if x[1] < y[1]: return -1
p.sort(c)
def namelist():
buffer = ""
for p in properties:
buffer += "calc_"+p[1]+", &\n"
buffer = buffer[:-4]
result = " namelist /properties/"+buffer
return result
def touch_all():
print "TOUCH",
for p in properties:
print "calc_"+p[1],
print ""
#file = open('../scripts/properties.py','w')
#print >>file,'properties = ',properties
#file.close()
def make_dims():
template = """
BEGIN_PROVIDER [ integer, size_%(p)s ]
implicit none
BEGIN_DOC
! Size of %(p)s
END_DOC
if (calc_%(p)s) then
size_%(p)s = %(d)s
else
size_%(p)s = 1
endif
END_PROVIDER
"""
for p in dims:
print template%{'p': p, 'd': dims[p]}

383
src/pseudo.irp.f Normal file
View File

@ -0,0 +1,383 @@
BEGIN_PROVIDER [ double precision, pseudo_v_k , (nucl_num,pseudo_klocmax) ]
implicit none
BEGIN_DOC
! V_k
END_DOC
call get_pseudo_pseudo_v_k(pseudo_v_k)
END_PROVIDER
BEGIN_PROVIDER [ double precision, pseudo_v_kl , (nucl_num,pseudo_kmax,0:pseudo_lmax) ]
implicit none
BEGIN_DOC
! V_Kl
END_DOC
call get_pseudo_pseudo_v_kl(pseudo_v_kl)
END_PROVIDER
BEGIN_PROVIDER [ integer, pseudo_kmax ]
implicit none
BEGIN_DOC
! kmax
END_DOC
call get_pseudo_pseudo_kmax(pseudo_kmax)
END_PROVIDER
BEGIN_PROVIDER [ double precision, pseudo_dz_kl , (nucl_num,pseudo_kmax,0:pseudo_lmax) ]
implicit none
BEGIN_DOC
! Exponents in the non-local part of the pseudo-potential
END_DOC
call get_pseudo_pseudo_dz_kl(pseudo_dz_kl)
END_PROVIDER
BEGIN_PROVIDER [ integer, pseudo_klocmax ]
implicit none
BEGIN_DOC
! klocmax
END_DOC
call get_pseudo_pseudo_klocmax(pseudo_klocmax)
END_PROVIDER
BEGIN_PROVIDER [ integer, pseudo_lmax ]
implicit none
BEGIN_DOC
! Max value of l in the pseudo-potential
END_DOC
call get_pseudo_pseudo_lmax(pseudo_lmax)
END_PROVIDER
BEGIN_PROVIDER [ integer, pseudo_grid_size ]
implicit none
BEGIN_DOC
! Size of the QMC grid (number of points)
END_DOC
call get_pseudo_pseudo_grid_size(pseudo_grid_size)
END_PROVIDER
BEGIN_PROVIDER [ double precision, pseudo_grid_rmax ]
implicit none
BEGIN_DOC
! Size of the QMC grid (max distance)
END_DOC
call get_pseudo_pseudo_grid_rmax(pseudo_grid_rmax)
END_PROVIDER
BEGIN_PROVIDER [ integer, pseudo_non_loc_dim ]
&BEGIN_PROVIDER [ integer, pseudo_non_loc_dim_8 ]
&BEGIN_PROVIDER [ integer, pseudo_non_loc_dim_count, (nucl_num) ]
implicit none
BEGIN_DOC
! Dimension of of pseudo_non_local arrays
END_DOC
pseudo_non_loc_dim = 0
pseudo_non_loc_dim_count = 0
integer :: k,l,m
do k=1,nucl_num
do l=0,pseudo_lmax
pseudo_non_loc_dim_count(k) += 2*l+1
enddo
enddo
pseudo_non_loc_dim = sum(pseudo_non_loc_dim_count)
integer, external :: mod_align
pseudo_non_loc_dim_8 = mod_align(pseudo_non_loc_dim)
END_PROVIDER
BEGIN_PROVIDER [ integer, pseudo_n_kl , (nucl_num,pseudo_kmax,0:pseudo_lmax) ]
implicit none
BEGIN_DOC
! n_kl
END_DOC
call get_pseudo_pseudo_n_kl(pseudo_n_kl)
END_PROVIDER
BEGIN_PROVIDER [ double precision, pseudo_dz_k , (nucl_num,pseudo_klocmax) ]
implicit none
BEGIN_DOC
! dz_k
END_DOC
call get_pseudo_pseudo_dz_k(pseudo_dz_k)
END_PROVIDER
BEGIN_PROVIDER [ integer, pseudo_n_k , (nucl_num,pseudo_klocmax) ]
implicit none
BEGIN_DOC
! n_k
END_DOC
call get_pseudo_pseudo_n_k(pseudo_n_k)
END_PROVIDER
BEGIN_PROVIDER [ logical, do_pseudo ]
implicit none
BEGIN_DOC
! Using pseudo potential integral of not
END_DOC
call get_pseudo_do_pseudo(do_pseudo)
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_pseudo_grid, (ao_num, -pseudo_lmax:pseudo_lmax, 0:pseudo_lmax, nucl_num, pseudo_grid_size) ]
implicit none
BEGIN_DOC
! Pseudopotential grid points
END_DOC
call get_pseudo_ao_pseudo_grid(ao_pseudo_grid)
END_PROVIDER
BEGIN_PROVIDER [ double precision, mo_pseudo_grid, (ao_num, -pseudo_lmax:pseudo_lmax, 0:pseudo_lmax, nucl_num, pseudo_grid_size) ]
implicit none
BEGIN_DOC
! Pseudopotential grid points
END_DOC
call get_pseudo_mo_pseudo_grid(mo_pseudo_grid)
END_PROVIDER
BEGIN_PROVIDER [ double precision, mo_pseudo_grid_scaled, (pseudo_non_loc_dim_8,ao_num,pseudo_grid_size) ]
implicit none
BEGIN_DOC
! Pseudopotential grid points
END_DOC
integer :: i,k,l,m,kk,n
double precision :: c
c = 1.d0/mo_scale
do n=1,pseudo_grid_size
do i=1,ao_num
kk = 0
do k=1,nucl_num
do l=0,pseudo_lmax
do m=-l,l
kk = kk+1
mo_pseudo_grid_scaled(kk,i,n) = c * mo_pseudo_grid(i,m,l,k,n)
enddo
enddo
enddo
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, v_pseudo_local, (elec_num) ]
implicit none
BEGIN_DOC
! Local component of the pseudo-potential
END_DOC
integer :: i,k,l
double precision :: r, rn, alpha
do l=1,elec_num
v_pseudo_local(l) = 0.d0
do k=1,pseudo_klocmax
do i=1,nucl_num
r = nucl_elec_dist(i,l)
alpha = pseudo_dz_k(i,k)*r*r
if (alpha > 20.d0) then
cycle
endif
select case (pseudo_n_k(i,k))
case (-2)
rn = nucl_elec_dist_inv(i,l)*nucl_elec_dist_inv(i,l)
case (-1)
rn = nucl_elec_dist_inv(i,l)
case (0)
rn = 1.d0
case (1)
rn = r
case (2)
rn = r*r
case default
rn = r**pseudo_n_k(i,k)
end select
v_pseudo_local(l) += pseudo_v_k(i,k) * exp(-alpha) * rn
enddo
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, v_pseudo_non_local, (pseudo_non_loc_dim_8,elec_num) ]
implicit none
BEGIN_DOC
! Non-Local component of the pseudo-potential
END_DOC
integer :: i,k,j,l
integer :: kk,m
double precision :: r, rn, r2, alpha
double precision :: tmp(0:pseudo_lmax)
v_pseudo_non_local = 0.d0
do j=1,elec_num
kk = 0
do i=1,nucl_num
r = nucl_elec_dist(i,j)
r2 = r*r
tmp(0:pseudo_lmax) = 0.d0
do k=1,pseudo_kmax
do l=0,pseudo_lmax
alpha = pseudo_dz_kl(i,k,l)*r2
if (alpha > 20.d0) then
cycle
endif
select case (pseudo_n_kl(i,k,l))
case (0)
rn = 1.d0
case (-1)
rn = nucl_elec_dist_inv(i,j)
case (1)
rn = r
case (2)
rn = r*r
case (-2)
rn = nucl_elec_dist_inv(i,j)*nucl_elec_dist_inv(i,j)
case default
rn = r**pseudo_n_kl(i,k,l)
end select
tmp(l) += pseudo_v_kl(i,k,l) * exp(-alpha) * rn
enddo
enddo
do l=0,pseudo_lmax
do m=-l,l
kk += 1
v_pseudo_non_local(kk,j) = tmp(l)
enddo
enddo
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, pseudo_mo_term, (mo_num,elec_num) ]
implicit none
BEGIN_DOC
! \sum < Ylm | MO > x Ylm(r) x V_nl(r)
END_DOC
integer :: ii,i,j,l,m,k,n,kk
double precision :: r, dr_inv
double precision :: tmp(pseudo_non_loc_dim_8,mo_num)
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: tmp
integer, save :: ifirst = 0
if (ifirst == 0) then
pseudo_mo_term = 0.d0
ifirst = 1
endif
PROVIDE pseudo_ylm present_mos mo_pseudo_grid_scaled pseudo_non_loc_dim_count
dr_inv = dble(pseudo_grid_size)/pseudo_grid_rmax
do j=1,elec_num
tmp = 0.d0
kk=0
do k=1,nucl_num
r = nucl_elec_dist(k,j)
n = 1 + int(0.5+r*dr_inv)
if (n<=pseudo_grid_size) then
do ii=1,num_present_mos
i = present_mos(ii)
!DIR$ LOOP COUNT(4)
do l=kk+1,kk+pseudo_non_loc_dim_count(k)
tmp(l,i) = mo_pseudo_grid_scaled (l,i,n) * pseudo_ylm(l,j)
enddo
enddo
endif
kk = kk+pseudo_non_loc_dim_count(k)
enddo
do ii=1,num_present_mos
i = present_mos(ii)
pseudo_mo_term(i,j) = 0.d0
do kk=1,pseudo_non_loc_dim
pseudo_mo_term(i,j) = pseudo_mo_term(i,j) + tmp(kk,i) * v_pseudo_non_local(kk,j)
enddo
enddo
enddo
END_PROVIDER
real function ylm(l,m,x,y,z,r_inv)
implicit none
include 'constants.F'
BEGIN_DOC
! Y(l,m) evaluated at x,y,z
END_DOC
integer, intent(in) :: l,m
real, intent(in) :: x,y,z,r_inv
ylm = 0.5 * one_over_sqpi
select case(l)
case(0)
continue
case(1)
select case (m)
case (-1)
ylm = ylm * sq3 * y * r_inv
case (0)
ylm = ylm * sq3 * z * r_inv
case (1)
ylm = ylm * sq3 * x * r_inv
end select
! case(2)
! select case (m)
! case(-2)
! ylm = ylm * sqrt(15.) * x * y * r_inv * r_inv
! case(-1)
! ylm = ylm * sqrt(15.) * y * z * r_inv * r_inv
! case(0)
! ylm = 0.5 * ylm * sqrt(15.) * (2.*z*z - x*x - y*y) * r_inv * r_inv
! case(1)
! ylm = ylm * sqrt(15.) * z * x * r_inv * r_inv
! case(2)
! ylm = 0.5 * ylm * sqrt(15.) * (x*x - y*y) * r_inv * r_inv
! end select
case default
stop 'problem in Ylm of pseudo'
end select
end
BEGIN_PROVIDER [ double precision, pseudo_ylm, (pseudo_non_loc_dim_8,elec_num) ]
implicit none
BEGIN_DOC
! Y(l,m) evaluated for every electron position centered on every nuclei
END_DOC
integer :: i,j,l,m,kk
real, external :: ylm
integer, save :: ifirst = 0
if (ifirst == 0) then
pseudo_ylm = 0.d0
ifirst = 1
endif
do j=1,elec_num
kk = 0
do i=1,nucl_num
do l=0,pseudo_lmax
do m=-l,l
kk = kk+1
pseudo_ylm(kk,j) = ylm(l,m, &
nucl_elec_dist_vec(1,i,j), &
nucl_elec_dist_vec(2,i,j), &
nucl_elec_dist_vec(3,i,j), &
nucl_elec_dist_inv(i,j))
enddo
enddo
enddo
enddo
END_PROVIDER

102
src/psi.irp.f Normal file
View File

@ -0,0 +1,102 @@
BEGIN_PROVIDER [ double precision, psi_value ]
implicit none
BEGIN_DOC
! Value of the wave function
END_DOC
psi_value = psidet_value*jast_value
if (psi_value == 0.d0) then
call abrt(irp_here,"Value of the wave function is 0.")
endif
END_PROVIDER
BEGIN_PROVIDER [ double precision, psi_value_inv ]
implicit none
BEGIN_DOC
! 1./psi_value
END_DOC
psi_value_inv = 1.d0/psi_value
END_PROVIDER
BEGIN_PROVIDER [ double precision, psi_value_inv2 ]
implicit none
BEGIN_DOC
! 1./(psi_value)**2
END_DOC
psi_value_inv2 = psi_value_inv*psi_value_inv
END_PROVIDER
BEGIN_PROVIDER [ double precision, psi_grad_x, (elec_num_8) ]
&BEGIN_PROVIDER [ double precision, psi_grad_y, (elec_num_8) ]
&BEGIN_PROVIDER [ double precision, psi_grad_z, (elec_num_8) ]
implicit none
BEGIN_DOC
! Gradients of the wave function
END_DOC
integer :: j
!DIR$ VECTOR ALIGNED
!DIR$ LOOP COUNT (100)
do j=1,elec_num
psi_grad_x(j) = psi_grad_psi_inv_x(j)*psi_value
psi_grad_y(j) = psi_grad_psi_inv_y(j)*psi_value
psi_grad_z(j) = psi_grad_psi_inv_z(j)*psi_value
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, psi_lapl, (elec_num_8) ]
implicit none
BEGIN_DOC
! Laplacian of the wave function
END_DOC
integer :: i, j
!DIR$ VECTOR ALIGNED
!DIR$ LOOP COUNT (100)
do j=1,elec_num
psi_lapl(j) = jast_value*(psidet_grad_lapl(4,j) + psidet_value*jast_lapl_jast_inv(j) + 2.d0*(&
psidet_grad_lapl(1,j)*jast_grad_jast_inv_x(j) + &
psidet_grad_lapl(2,j)*jast_grad_jast_inv_y(j) + &
psidet_grad_lapl(3,j)*jast_grad_jast_inv_z(j) ))
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, psi_grad_psi_inv_x, (elec_num_8) ]
&BEGIN_PROVIDER [ double precision, psi_grad_psi_inv_y, (elec_num_8) ]
&BEGIN_PROVIDER [ double precision, psi_grad_psi_inv_z, (elec_num_8) ]
implicit none
BEGIN_DOC
! grad(psi)/psi
END_DOC
integer :: j
!DIR$ VECTOR ALIGNED
!DIR$ LOOP COUNT (100)
do j=1,elec_num
psi_grad_psi_inv_x(j) = psidet_grad_lapl(1,j)*psidet_inv + jast_grad_jast_inv_x(j)
psi_grad_psi_inv_y(j) = psidet_grad_lapl(2,j)*psidet_inv + jast_grad_jast_inv_y(j)
psi_grad_psi_inv_z(j) = psidet_grad_lapl(3,j)*psidet_inv + jast_grad_jast_inv_z(j)
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, psi_lapl_psi_inv, (elec_num_8) ]
implicit none
BEGIN_DOC
! (Laplacian psi) / psi
END_DOC
integer :: i, j
!DIR$ VECTOR ALIGNED
!DIR$ LOOP COUNT (100)
do j=1,elec_num
psi_lapl_psi_inv(j) = psi_lapl(j)*psi_value_inv
enddo
END_PROVIDER

350
src/simulation.irp.f Normal file
View File

@ -0,0 +1,350 @@
BEGIN_PROVIDER [ logical, is_worker ]
implicit none
BEGIN_DOC
! True if the process is a worker process
END_DOC
is_worker = .False.
END_PROVIDER
BEGIN_PROVIDER [ integer, walk_num_tot ]
implicit none
BEGIN_DOC
! Total number of walkers
END_DOC
walk_num_tot = 1000
call get_electrons_elec_walk_num_tot(walk_num_tot)
walk_num_tot = max(walk_num,walk_num_tot)
call iinfo(irp_here,'walk_num', walk_num_tot)
if (walk_num_tot <= 0) then
call abrt(irp_here,'Total number of walkers should be > 0')
endif
END_PROVIDER
BEGIN_PROVIDER [ integer, walk_num ]
&BEGIN_PROVIDER [ integer, walk_num_8 ]
implicit none
BEGIN_DOC
! Number of walkers
END_DOC
walk_num = 100
call get_electrons_elec_walk_num(walk_num)
call iinfo(irp_here,'walk_num', walk_num)
if (walk_num <= 0) then
call abrt(irp_here,'Number of walkers should be > 0')
endif
integer :: mod_align
walk_num_8 = mod_align(walk_num)
END_PROVIDER
BEGIN_PROVIDER [ logical, do_equilibration ]
implicit none
BEGIN_DOC
! Equilibrate walkers
END_DOC
do_equilibration = .True.
if (.not.do_prepare) then
call get_simulation_equilibration(do_equilibration)
endif
call iinfo(irp_here,'equilibration', do_equilibration)
END_PROVIDER
BEGIN_PROVIDER [ logical, do_prepare ]
implicit none
BEGIN_DOC
! If true, prepare new walkers
END_DOC
do_prepare = .False.
END_PROVIDER
BEGIN_PROVIDER [ double precision, block_time ]
implicit none
BEGIN_DOC
! Wall time requested to realize one block
END_DOC
block_time = 30.d0
integer :: block_time_int
call get_simulation_block_time(block_time_int)
if (block_time<= 1) then
call abrt(irp_here,'Block time should be > 1s')
endif
double precision, external :: qmc_ranf
block_time = dble(block_time_int) + qmc_ranf()
call dinfo(irp_here,'block_time',block_time)
END_PROVIDER
BEGIN_PROVIDER [ integer, stop_time ]
implicit none
BEGIN_DOC
! Termination condition of the run
END_DOC
stop_time = 3600*24
call get_simulation_stop_time(stop_time)
call iinfo(irp_here,'stop_time',stop_time)
if (stop_time<= 1) then
call abrt(irp_here,'Stop time should be > 1s')
endif
END_PROVIDER
BEGIN_PROVIDER [ real, time_step ]
&BEGIN_PROVIDER [ real, time_step_inv ]
&BEGIN_PROVIDER [ double precision, dtime_step ]
implicit none
BEGIN_DOC
! time_step : The time step of the random walk
END_DOC
time_step = 0.0
call get_simulation_time_step(time_step)
call rinfo(irp_here,'time_step',time_step)
if (time_step <= 0.) then
call abrt(irp_here,'Time step should be > 0')
endif
dtime_step = dble(time_step)
time_step_inv = 1./time_step
END_PROVIDER
BEGIN_PROVIDER [ double precision, time_step_sq ]
&BEGIN_PROVIDER [ double precision, time_step_exp ]
&BEGIN_PROVIDER [ double precision, time_step_exp_sq ]
&BEGIN_PROVIDER [ double precision, time_step_exp_sq_sq ]
implicit none
BEGIN_DOC
!
! time_step_sq : sqrt(time_step)
!
! time_step_exp = exp(-time_step)
!
! time_step_exp_sq = sqrt(time_step_exp)
!
! time_step_exp_sq_sq = sqrt(time_step_exp_sq)
END_DOC
time_step_sq = sqrt(dble(time_step))
time_step_exp = exp(-dble(time_step))
time_step_exp_sq = sqrt(time_step_exp)
time_step_exp_sq_sq = sqrt(time_step_exp_sq)
END_PROVIDER
BEGIN_PROVIDER [ integer, qmc_method ]
implicit none
include 'types.F'
BEGIN_DOC
! qmc_method : Calculation method. Can be t_VMC, t_DMC
END_DOC
character*(32) :: method
method = types(t_VMC)
call get_simulation_method(method)
if (method == types(t_VMC)) then
qmc_method = t_VMC
else if (method == types(t_DMC)) then
qmc_method = t_DMC
else
call abrt(irp_here, 'Method should be ( VMC | DMC )')
endif
call cinfo(irp_here,'qmc_method',trim(method))
END_PROVIDER
BEGIN_PROVIDER [ real, events_num ]
BEGIN_DOC
! Number of Monte Carlo events to average
END_DOC
events_num = real(walk_num)*real(step_num)
END_PROVIDER
BEGIN_PROVIDER [ integer, walk_i ]
BEGIN_DOC
! Current walker
END_DOC
walk_i = 1
END_PROVIDER
BEGIN_PROVIDER [ real, accepted_num ]
&BEGIN_PROVIDER [ real, rejected_num ]
BEGIN_DOC
! Number of accepted steps
! Number of rejected steps
END_DOC
accepted_num= 0.
rejected_num= 0.
END_PROVIDER
BEGIN_PROVIDER [ logical, save_data ]
implicit none
BEGIN_DOC
! If true, the updated simulation data is saved for restart.
END_DOC
save_data = .True.
call get_simulation_save_data(save_data)
call linfo(irp_here,'save_data',save_data)
END_PROVIDER
real function accep_rate()
if ((accepted_num+rejected_num) > 0.) then
accep_rate = accepted_num/(accepted_num+rejected_num)
else
accep_rate = 0.
endif
end
logical function first_step()
first_step = (accepted_num+rejected_num < 1.)
end
subroutine accep_reset
FREE accepted_num
FREE rejected_num
end
BEGIN_PROVIDER [ integer, print_level ]
BEGIN_DOC
! Level of verbosity for standard output printing
END_DOC
print_level = 1
call get_simulation_print_level(print_level)
END_PROVIDER
BEGIN_PROVIDER [ character*(64), hostname]
implicit none
BEGIN_DOC
! Name of the current host
END_DOC
call HOSTNM(hostname)
END_PROVIDER
BEGIN_PROVIDER [ logical, do_nucl_fitcusp ]
implicit none
BEGIN_DOC
! If true, do the fit of the electron-nucleus cusp
END_DOC
do_nucl_fitcusp = .True.
call get_simulation_do_nucl_fitcusp(do_nucl_fitcusp)
call linfo(irp_here,'do_nucl_fitcusp',do_nucl_fitcusp)
END_PROVIDER
BEGIN_PROVIDER [ integer, vmc_algo ]
implicit none
include 'types.F'
BEGIN_DOC
! Type of VMC algorithm: Brownian, MTM or Langevin
END_DOC
character*(32) :: Sampling
Sampling = types(t_Langevin)
call get_simulation_sampling(Sampling)
vmc_algo = 0
if (Sampling == types(t_Brownian)) then
vmc_algo = t_Brownian
else if (Sampling == types(t_Langevin)) then
vmc_algo = t_Langevin
if (qmc_method == t_DMC) then
vmc_algo = t_Brownian
endif
else if (Sampling == types(t_MTM)) then
vmc_algo = t_MTM
else
call abrt(irp_here,'Sampling should be (Brownian|Langevin|MTM|Read)')
endif
call cinfo(irp_here,'Sampling',Sampling)
ASSERT (vmc_algo > 0)
END_PROVIDER
BEGIN_PROVIDER [ character*(512), ezfio_filename ]
implicit none
BEGIN_DOC
! Name of the ezfio file.
! Defined in init_ezfio_filename
END_DOC
integer :: command_argument_count
if (command_argument_count() == 0) then
ezfio_filename = 'NOT_SET'
call ezfio_set_file(ezfio_filename)
else
call get_command_argument(1,ezfio_filename)
if (.not.is_worker) then
call ezfio_set_file(ezfio_filename)
endif
endif
END_PROVIDER
BEGIN_PROVIDER [ character*(128), http_server ]
implicit none
BEGIN_DOC
! Address of the data server
END_DOC
integer :: command_argument_count
if (command_argument_count() > 1) then
call get_command_argument(2,http_server)
else
call get_simulation_http_server(http_server)
endif
END_PROVIDER
subroutine read_do_run(do_run)
implicit none
integer, intent(out) :: do_run
BEGIN_DOC
! Read the do_run variable from the ezfio directory
END_DOC
include 'types.F'
do_run = t_Stopping
call get_simulation_do_run(do_run)
end
BEGIN_PROVIDER [ character*(32), md5_key ]
implicit none
BEGIN_DOC
! Digest of the input
END_DOC
md5_key = ''
call get_simulation_md5_key(md5_key)
if (md5_key == '') then
call abrt(irp_here,'MD5 key of input is absent')
endif
END_PROVIDER

34
src/types.F Normal file
View File

@ -0,0 +1,34 @@
integer, parameter :: t_Brownian = 3
integer, parameter :: t_Langevin = 4
integer, parameter :: t_MTM = 5
integer, parameter :: t_Read = 6
integer, parameter :: t_VMC = 7
integer, parameter :: t_DMC = 8
integer, parameter :: t_Simple = 11
integer, parameter :: t_None = 12
integer, parameter :: t_Core = 14
integer, parameter :: t_Stopped = 0
integer, parameter :: t_Queued = 1
integer, parameter :: t_Running = 2
integer, parameter :: t_Stopping = 3
character*(32) :: types(15) = &
(/ '', &
'', &
'Brownian', &
'Langevin', &
'', &
'', &
'VMC ', &
'DMC ', &
' ', &
'', &
'Simple ', &
'None ', &
' ', &
'Core ', &
' '/)

370
src/wf.irp.f Normal file
View File

@ -0,0 +1,370 @@
BEGIN_PROVIDER [ integer, i_state ]
implicit none
BEGIN_DOC
! Current state
END_DOC
i_state = 1
END_PROVIDER
BEGIN_PROVIDER [ integer, N_int ]
implicit none
BEGIN_DOC
! Number of 64-bit integers needed to represent determinants as binary strings
END_DOC
call get_spindeterminants_n_int(N_int)
END_PROVIDER
BEGIN_PROVIDER [ integer, bit_kind ]
implicit none
BEGIN_DOC
! Number of octets per integer storing determinants
END_DOC
call get_spindeterminants_bit_kind(bit_kind)
ASSERT (bit_kind == 8)
END_PROVIDER
BEGIN_PROVIDER [ integer, N_states ]
implicit none
BEGIN_DOC
! Number of states in EZFIO file
END_DOC
call get_spindeterminants_n_states(N_states)
END_PROVIDER
BEGIN_PROVIDER [ integer, det_num_input ]
implicit none
BEGIN_DOC
! Number of Det_a x Det_b products in input file
END_DOC
call get_spindeterminants_n_det(det_num_input)
END_PROVIDER
BEGIN_PROVIDER [ double precision, det_alpha_norm, (det_alpha_num) ]
&BEGIN_PROVIDER [ double precision, det_beta_norm, (det_beta_num) ]
implicit none
BEGIN_DOC
! Norm of the alpha and beta spin determinants in the wave function:
!
! ||Da||_i \sum_j C_{ij}**2
END_DOC
integer :: i,j,k
double precision :: f
det_alpha_norm = 0.d0
det_beta_norm = 0.d0
do k=1,det_num
i = det_coef_matrix_rows(k)
j = det_coef_matrix_columns(k)
f = det_coef_matrix_values(k)*det_coef_matrix_values(k)
det_alpha_norm(i) += f
det_beta_norm(j) += f
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, det_coef_matrix_values, (det_num_input) ]
&BEGIN_PROVIDER [ integer, det_coef_matrix_rows, (det_num_input) ]
&BEGIN_PROVIDER [ integer, det_coef_matrix_columns, (det_num_input) ]
implicit none
BEGIN_DOC
! det_coef_matrix in sparse storage (Coordinate format for sparse BLAS)
END_DOC
double precision, allocatable :: buffer(:,:)
allocate (buffer(det_num_input,N_states))
call get_spindeterminants_psi_coef_matrix_rows(det_coef_matrix_rows)
call get_spindeterminants_psi_coef_matrix_columns(det_coef_matrix_columns)
call get_spindeterminants_psi_coef_matrix_values(buffer)
det_coef_matrix_values(:) = buffer(:,i_state)
deallocate(buffer)
END_PROVIDER
BEGIN_PROVIDER [ integer, det_num ]
implicit none
BEGIN_DOC
! Number of Det_a x Det_b products. The determinant basis set is reduced with
! the CI threshold
END_DOC
integer :: i,j,k,l
double precision :: f
double precision :: d_alpha(det_alpha_num), d_beta (det_beta_num)
integer :: i_alpha(det_alpha_num), i_beta(det_beta_num)
integer :: iorder(max(det_alpha_num,det_beta_num))
double precision :: t, norm
t = ci_threshold
! Compute the norm of the alpha and beta determinants
d_alpha = 0.d0
d_beta = 0.d0
do k=1,det_num_input
i = det_coef_matrix_rows(k)
j = det_coef_matrix_columns(k)
f = det_coef_matrix_values(k)*det_coef_matrix_values(k)
d_alpha(i) += f
d_beta (j) += f
enddo
t = min(t, maxval(d_alpha))
t = min(t, maxval(d_beta))
! Reorder alpha determinants
do i=1,det_alpha_num
iorder(i) = i
if (d_alpha(i) < t) then
i_alpha(i) = det_alpha_num+i
else
i_alpha(i) = i
endif
enddo
call isort(i_alpha,iorder,det_alpha_num)
i=det_alpha_num
do while (i > 0)
if (i_alpha(i) <= det_alpha_num) then
det_alpha_num = i
exit
else
i = i-1
endif
enddo
do i=1,det_alpha_num
psi_det_alpha(:,i) = psi_det_alpha(:,iorder(i))
i_alpha(iorder(i)) = i
enddo
! Reorder beta determinants
do i=1,det_beta_num
iorder(i) = i
if (d_beta(i) < t) then
i_beta(i) = det_beta_num+i
else
i_beta(i) = i
endif
enddo
call isort(i_beta,iorder,det_beta_num)
i=det_beta_num
do while (i > 0)
if (i_beta(i) <= det_beta_num) then
det_beta_num = i
exit
else
i = i-1
endif
enddo
do i=1,det_beta_num
psi_det_beta(:,i) = psi_det_beta(:,iorder(i))
i_beta(iorder(i)) = i
enddo
! Apply the threshold to the wave function
l = 1
norm = 0.d0
do k=1,det_num_input
i = det_coef_matrix_rows(k)
j = det_coef_matrix_columns(k)
det_coef_matrix_rows(l) = i_alpha(i)
det_coef_matrix_columns(l) = i_beta(j)
det_coef_matrix_values(l) = det_coef_matrix_values(k)
if ( (d_alpha(i) >= t).and.(d_beta(j) >= t) ) then
l = l+1
norm += det_coef_matrix_values(k)*det_coef_matrix_values(k)
endif
enddo
det_num = l-1
norm = 1.d0/dsqrt(norm)
do k=1,det_num
det_coef_matrix_values(k) *= norm
enddo
SOFT_TOUCH det_alpha_num det_beta_num det_coef_matrix_values det_coef_matrix_rows det_coef_matrix_columns psi_det_beta psi_det_alpha
END_PROVIDER
BEGIN_PROVIDER [ integer, det_alpha_num ]
&BEGIN_PROVIDER [ integer, det_beta_num ]
implicit none
BEGIN_DOC
! Number of alpha and beta determinants
END_DOC
call get_spindeterminants_n_det_alpha(det_alpha_num)
call get_spindeterminants_n_det_beta(det_beta_num)
END_PROVIDER
BEGIN_PROVIDER [ integer, det_alpha_num_8 ]
&BEGIN_PROVIDER [ integer, det_beta_num_8 ]
implicit none
BEGIN_DOC
! Number of alpha and beta determinants
END_DOC
integer :: mod_align
det_alpha_num_8 = max(4,mod_align(det_alpha_num)) !
det_beta_num_8 = max(4,mod_align(det_beta_num)) ! Used in 4x unrolling
END_PROVIDER
BEGIN_PROVIDER [ double precision, ci_threshold ]
implicit none
BEGIN_DOC
! Threshold on absolute value of the CI coefficients of the wave functioE
END_DOC
ci_threshold = 0.d0
call get_simulation_ci_threshold(ci_threshold)
call rinfo(irp_here,'ci_threshold',ci_threshold)
END_PROVIDER
BEGIN_PROVIDER [ integer*8, psi_det_alpha, (N_int,det_alpha_num) ]
implicit none
BEGIN_DOC
! Alpha determinants
END_DOC
call get_spindeterminants_psi_det_alpha(psi_det_alpha)
END_PROVIDER
BEGIN_PROVIDER [ integer*8, psi_det_beta, (N_int,det_beta_num) ]
implicit none
BEGIN_DOC
! Beta determinants
END_DOC
call get_spindeterminants_psi_det_beta(psi_det_beta)
END_PROVIDER
BEGIN_PROVIDER [ integer, present_mos, (mo_tot_num) ]
&BEGIN_PROVIDER [ integer, num_present_mos ]
&BEGIN_PROVIDER [ integer, num_present_mos_8 ]
&BEGIN_PROVIDER [ integer, mo_closed_num ]
implicit none
BEGIN_DOC
! List of used MOs to build the wf in the CI expansion
END_DOC
integer*8 :: tmp_det(N_int)
integer :: i,k
integer, external :: mod_align
PROVIDE det_num
num_present_mos = mo_tot_num
do i=1,mo_tot_num
present_mos(i) = i
enddo
!---
present_mos = 0
tmp_det = 0_8
do i=1,det_alpha_num
do k=1,N_int
tmp_det(k) = ior(tmp_det(k),psi_det_alpha(k,i))
enddo
enddo
do i=1,det_beta_num
do k=1,N_int
tmp_det(k) = ior(tmp_det(k),psi_det_beta(k,i))
enddo
enddo
call bitstring_to_list(tmp_det,present_mos,num_present_mos,N_int)
!---
num_present_mos_8 = mod_align(num_present_mos)
integer :: list(mo_tot_num), n
logical :: good
list = present_mos
mo_closed_num = elec_beta_num
do n=1,elec_beta_num
call list_to_bitstring(tmp_det,present_mos,n,N_int)
do k=1,N_int
if (tmp_det(k) == 0_8) then
exit
endif
good = .True.
do i=1,det_alpha_num
if (iand(tmp_det(k),psi_det_alpha(k,i)) /= tmp_det(k)) then
good = .False.
exit
endif
enddo
if (good) then
do i=1,det_beta_num
if (iand(tmp_det(k),psi_det_beta(k,i)) /= tmp_det(k)) then
good = .False.
exit
endif
enddo
endif
if (.not.good) then
exit
endif
enddo
if (.not.good) then
mo_closed_num = n-1
exit
endif
enddo
END_PROVIDER
subroutine list_to_bitstring( string, list, n_elements, Nint)
implicit none
BEGIN_DOC
! Returns the physical string "string(N_int,2)" from the array of
! occupations "list(N_int*64,2)
END_DOC
integer, intent(in) :: Nint
integer*8, intent(out) :: string(Nint)
integer, intent(in) :: list(Nint*64)
integer, intent(in) :: n_elements
integer :: i, j
integer :: ipos, iint
!
! <== ipos ==>
! |
! v
!string :|------------------------|-------------------------|------------------------|
! <==== 64 ====> <==== 64 ====> <==== 64 ====>
! { iint } { iint } { iint }
!
string = 0_8
do i=1,n_elements
iint = ishft(list(i)-1,-6) + 1
ipos = list(i)-ishft((iint-1),6)-1
string(iint) = ibset( string(iint), ipos )
enddo
end
BEGIN_PROVIDER [ integer, det_alpha_order, (det_alpha_num) ]
implicit none
BEGIN_DOC
! Order in which to compute the alhpa determinants
END_DOC
integer :: i
! double precision :: tmp(det_alpha_num)
do i=1,det_alpha_num
det_alpha_order(i) = i
enddo
END_PROVIDER
BEGIN_PROVIDER [ integer, det_beta_order, (det_beta_num) ]
implicit none
BEGIN_DOC
! Order in which to compute the beta determinants
END_DOC
integer :: i
do i=1,det_beta_num
det_beta_order(i) = i
enddo
END_PROVIDER