LCPQ/qmcchem
10
1
Fork 1
mirror of https://gitlab.com/scemama/qmcchem.git synced 2024-12-21 11:53:30 +01:00
Code Releases Activity

Merge branch 'master' of gitlab.com:scemama/qmcchem

Browse Source
This commit is contained in:
Anthony Scemama 2021-08-26 11:41:45 +02:00
commit 6617cf09a1
11 changed files with 409 additions and 10 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
bin/jast_opt.py
View File

@ -50,6 +50,10 @@ def main():
def get_params_pen():
d = ezfio.jastrow_jast_pen
print(atom_map)
for m in atom_map:
print(m[0])
print (d[m[0]])
return np.array([d[m[0]] for m in atom_map])
@ -137,8 +141,8 @@ def main():
variance, v_err = get_variance()
if e is None or variance is None:
continue
energy = e + variance
err = sqrt(e_err*e_err+v_err*v_err)
energy = e #+ variance
err = e_err #sqrt(e_err*e_err+v_err*v_err)
print(" %f %f %f %f %f %f"%(e, e_err, variance, v_err, energy, err))
if (energy-2.*err) > memo_energy['fmin']+thresh:
local_thresh = 10.*thresh

1
ezfio_config/qmc.config
View File

@ -83,6 +83,7 @@ jastrow
jast_a_up_dn real
jast_b_up_up real
jast_b_up_dn real
mu_erf real
jast_pen real (nuclei_nucl_num)
jast_eeN_e_a real (nuclei_nucl_num)
jast_eeN_e_b real (nuclei_nucl_num)

10
ocaml/Input.ml
View File

@ -775,7 +775,7 @@ end
module Jastrow_type : sig
type t = None | Core | Simple
type t = None | Core | Simple | Mu
val doc : string
val read : unit -> t
val write : t -> unit
@ -784,20 +784,22 @@ module Jastrow_type : sig
end = struct
type t = None | Core | Simple
let doc = "Type of Jastrow factor [ None | Core | Simple ]"
type t = None | Core | Simple | Mu
let doc = "Type of Jastrow factor [ None | Core | Simple | Mu ]"
let of_string s =
match String.capitalize_ascii (String.trim s) with
| "Core" -> Core
| "Simple" -> Simple
| "None" -> None
| _ -> failwith "Jastrow type should be [ None | Core | Simple ]"
| "Mu" -> Mu
| _ -> failwith "Jastrow type should be [ None | Core | Simple | Mu ]"
let to_string = function
| Core -> "Core"
| Simple -> "Simple"
| Mu -> "Mu"
| None -> "None"

3
src/JASTROW/jastrow_full.irp.f
View File

@ -26,6 +26,7 @@ BEGIN_TEMPLATE
SUBST [X]
Simple ;;
Core ;;
Mu ;;
END_TEMPLATE
if (ifirst == 0) then
dshift = argexpo
@ -83,6 +84,7 @@ BEGIN_TEMPLATE
SUBST [ X ]
Simple ;;
Core ;;
Mu ;;
END_TEMPLATE
!DIR$ VECTOR ALIGNED
!DIR$ LOOP COUNT (200)
@ -129,6 +131,7 @@ BEGIN_TEMPLATE
SUBST [X]
Simple ;;
Core ;;
Mu ;;
END_TEMPLATE
!DIR$ VECTOR ALIGNED

128
src/JASTROW/jastrow_mu.irp.f Normal file
View File

@ -0,0 +1,128 @@
! Mu Jastrow
! --------------
! See Giner JCP 2021
BEGIN_PROVIDER [ double precision , jast_elec_Mu_value, (elec_num_8) ]
implicit none
BEGIN_DOC
! J(i) = \sum_j a.rij/(1+b^2.rij) - \sum_A (a.riA/(1+a.riA))^2
! Eq (11)
END_DOC
integer :: i,j
double precision :: a, b, rij, tmp
include '../constants.F'
double precision :: mu
mu = mu_erf
do i=1,elec_num
jast_elec_Mu_value(i) = 0.d0
enddo
do j=1,elec_num
!DIR$ LOOP COUNT (50)
do i=1,elec_num
if(j==i)cycle
rij = elec_dist(i,j)
tmp = 0.5d0 * rij * (1.d0 - derf(mu*rij)) - 0.5d0/(dsqpi*mu) * dexp(-mu*mu*rij*rij)
jast_elec_Mu_value(i) += tmp
enddo
enddo
jast_elec_Mu_value = jast_elec_Mu_value * 0.5d0 ! symmetrization
END_PROVIDER
BEGIN_PROVIDER [ double precision , jast_elec_Mu_grad_x, (elec_num_8) ]
&BEGIN_PROVIDER [ double precision , jast_elec_Mu_grad_y, (elec_num_8) ]
&BEGIN_PROVIDER [ double precision , jast_elec_Mu_grad_z, (elec_num_8) ]
implicit none
BEGIN_DOC
! Gradient of the Jastrow factor
! Eq (A1)
END_DOC
integer :: i,j
double precision :: a, b, rij, tmp, x, y, z
include '../constants.F'
double precision :: mu
mu = mu_erf
do i=1,elec_num
jast_elec_Mu_grad_x(i) = 0.d0
jast_elec_Mu_grad_y(i) = 0.d0
jast_elec_Mu_grad_z(i) = 0.d0
!DIR$ LOOP COUNT (100)
enddo
! (grad of J(r12) with respect to xi, yi, zi)
do i = 1, elec_num
do j = 1, elec_num
if(i==j)cycle
rij = elec_dist(j,i)
jast_elec_Mu_grad_x(i) += 0.5d0 * ( 1.d0 - derf(mu * rij) ) * elec_dist_inv(j,i) * (-1.d0) * elec_dist_vec_x(j,i)
jast_elec_Mu_grad_y(i) += 0.5d0 * ( 1.d0 - derf(mu * rij) ) * elec_dist_inv(j,i) * (-1.d0) * elec_dist_vec_y(j,i)
jast_elec_Mu_grad_z(i) += 0.5d0 * ( 1.d0 - derf(mu * rij) ) * elec_dist_inv(j,i) * (-1.d0) * elec_dist_vec_z(j,i)
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision , jast_elec_Mu_lapl, (elec_num_8) ]
implicit none
BEGIN_DOC
! Laplacian of the Jastrow factor
! Eq (A10)
END_DOC
integer :: i,j
double precision :: a, b, rij, tmp, x, y, z
include '../constants.F'
double precision :: mu, x_ij, y_ij, z_ij, rij_inv
mu = mu_erf
do i=1,elec_num
jast_elec_Mu_lapl(i) = 0.d0
enddo
do i=1, elec_num
do j=1, elec_num
if(j==i)cycle
rij = elec_dist(j,i)
rij_inv = elec_dist_inv(j,i)
x_ij = elec_dist_vec_x(j,i)
y_ij = elec_dist_vec_y(j,i)
z_ij = elec_dist_vec_z(j,i)
jast_elec_Mu_lapl(i) += (1.d0 - derf(mu*rij))*elec_dist_inv(j,i) - mu/dsqpi * dexp(-mu*mu*rij*rij)
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, mu_erf ]
implicit none
mu_erf = 0.5d0
END_PROVIDER
BEGIN_PROVIDER [double precision, grad_j_mu_x,(elec_num, elec_num)]
&BEGIN_PROVIDER [double precision, grad_j_mu_y,(elec_num, elec_num)]
&BEGIN_PROVIDER [double precision, grad_j_mu_z,(elec_num, elec_num)]
implicit none
BEGIN_DOC
! Needed for 3-body terms
END_DOC
integer :: i,j
double precision :: rij, mu,scal
mu = mu_erf
grad_j_mu_x = 0.d0
grad_j_mu_y = 0.d0
grad_j_mu_z = 0.d0
do j = 1, elec_num
do i = 1, elec_num
if(i==j)cycle
rij = elec_dist(i,j)
scal = 0.5d0 * ( 1.d0 - derf(mu * rij) ) * elec_dist_inv(i,j)
grad_j_mu_x(i,j) = (elec_coord_transp(1,i) - elec_coord_transp(1,j)) * scal
grad_j_mu_y(i,j) = (elec_coord_transp(2,i) - elec_coord_transp(2,j)) * scal
grad_j_mu_z(i,j) = (elec_coord_transp(3,i) - elec_coord_transp(3,j)) * scal
enddo
enddo
END_PROVIDER

4
src/JASTROW/jastrow_param.irp.f
View File

@ -27,8 +27,10 @@ BEGIN_PROVIDER [ integer, jast_type ]
jast_type = t_None
else if (buffer == types(t_Core)) then
jast_type = t_Core
else if (buffer == types(t_Mu)) then
jast_type = t_Mu
else
call abrt(irp_here,'Jastrow type should be (None|Simple|Core)')
call abrt(irp_here,'Jastrow type should be (None|Simple|Core|Mu)')
endif
call cinfo(irp_here,'jast_type',buffer)

25
src/MAIN/vmc_test.irp.f Normal file
View File

@ -0,0 +1,25 @@
program vmc_test
real :: t1,t2
print *, 'Ndet=',det_num
print *, 'Ndet alpha beta =',det_alpha_num, det_beta_num
if (do_prepare) then
stop 'No walkers'
endif
print *, 'E_loc = ', E_loc
call step2
call ezfio_finish
end
subroutine step2
implicit none
real :: accep_rate
print *, '---'
print *, '<E_loc> = ', E_loc_block_walk
print *, '<E_loc_2> = ', E_loc_2_block_walk
print *, 'w = ', block_weight
print *, 'Accept', accep_rate()
print *, ''
print *, ci_dress_vmc_block_walk
end

72
src/PROPERTIES/properties_ci.irp.f
View File

@ -1,3 +1,15 @@
BEGIN_PROVIDER [ double precision, psi_norm ]
implicit none
BEGIN_DOC
! <1/J^2>
END_DOC
psi_norm = jast_value_inv*jast_value_inv
psi_norm_min = min(psi_norm_min,psi_norm)
psi_norm_max = max(psi_norm_max,psi_norm)
SOFT_TOUCH psi_norm_min psi_norm_max
END_PROVIDER
BEGIN_PROVIDER [ double precision, ci_overlap_psidet, (size_ci_overlap_psidet) ]
implicit none
BEGIN_DOC
@ -29,7 +41,7 @@ BEGIN_PROVIDER [ double precision, ci_h_psidet, (size_ci_h_psidet) ]
END_DOC
integer :: i, j, k, l
double precision :: T, tmp
double precision :: T
do k=1,det_num
i = det_coef_matrix_rows(k)
@ -42,7 +54,7 @@ BEGIN_PROVIDER [ double precision, ci_h_psidet, (size_ci_h_psidet) ]
T += det_beta_grad_lapl (4,l,j)*det_alpha_value(i)
enddo
ci_h_psidet(k) = -0.5d0*T + E_pot * det_alpha_value(i)*det_beta_value (j)
ci_h_psidet(k) *= psidet_inv
ci_h_psidet(k) *= psi_value_inv * jast_value_inv
enddo
ci_h_psidet_min = min(ci_h_psidet_min,minval(ci_h_psidet))
@ -217,3 +229,59 @@ BEGIN_PROVIDER [ double precision, ci_h_matrix_diag, (size_ci_h_matrix_diag) ]
END_PROVIDER
BEGIN_PROVIDER [ double precision, ci_h_transcor_psi, (size_ci_h_transcor_psi) ]
implicit none
BEGIN_DOC
! < det(i) e^{-J} |H| Psi >
!
! Dimensions : det_num
END_DOC
integer :: i, j, k
do k=1,det_num
i = det_coef_matrix_rows(k)
j = det_coef_matrix_columns(k)
ci_h_transcor_psi(k) = E_loc * jast_value_inv * &
det_alpha_value(i)*det_beta_value(j) * psi_value_inv
enddo
ci_h_transcor_psi_min = min(ci_h_transcor_psi_min,minval(ci_h_transcor_psi))
ci_h_transcor_psi_max = max(ci_h_transcor_psi_max,maxval(ci_h_transcor_psi))
SOFT_TOUCH ci_h_transcor_psi_min ci_h_transcor_psi_max
END_PROVIDER
BEGIN_PROVIDER [ double precision, ci_dress, (size_ci_dress) ]
implicit none
BEGIN_DOC
! < det(i) e^{-J} |H| Psi >
!
! Dimensions : det_num
END_DOC
integer :: i, j, k, l
double precision :: T, h_psidet
do k=1,det_num
i = det_coef_matrix_rows(k)
j = det_coef_matrix_columns(k)
T = 0.d0
do l=1,elec_alpha_num
T += det_alpha_grad_lapl(4,l,i)*det_beta_value (j)
enddo
do l=1,elec_beta_num
T += det_beta_grad_lapl (4,l,j)*det_alpha_value(i)
enddo
h_psidet = -0.5d0*T + E_pot * det_alpha_value(i)*det_beta_value(j)
h_psidet *= psi_value_inv * jast_value_inv
ci_dress(k) = E_loc * det_alpha_value(i)*det_beta_value(j) * psi_value_inv * jast_value_inv &
- h_psidet
enddo
ci_dress_min = min(ci_dress_min,minval(ci_dress))
ci_dress_max = max(ci_dress_max,maxval(ci_dress))
SOFT_TOUCH ci_dress_min ci_dress_max
END_PROVIDER

164
src/PROPERTIES/properties_mu.irp.f Normal file
View File

@ -0,0 +1,164 @@
BEGIN_PROVIDER [ double precision, Energy_mu ]
implicit none
BEGIN_DOC
! E mu
! Eq.(26-30)
END_DOC
integer :: i
energy_mu = E_nucl
!DIR$ VECTOR ALIGNED
!DIR$ LOOP COUNT(200)
do i=1,elec_num
energy_mu += E_kin_elec(i)
enddo
energy_mu += Eff_pot_mu + eff_pot_deriv_mu + E_nucl_elec - three_body_mu
energy_mu_min = min(energy_mu_min,energy_mu)
energy_mu_max = max(energy_mu_max,energy_mu)
SOFT_TOUCH energy_mu_min energy_mu_max
END_PROVIDER
BEGIN_PROVIDER [double precision, E_nucl_elec]
implicit none
!TODO
integer :: i,j
E_nucl_elec = 0.d0
do i = 1, elec_num
! E_nucl_elec += E_pot_elec_one(i) + E_pot_elec_two(i)
E_nucl_elec += E_pot_elec_one(i)
enddo
E_nucl_elec_min = min(E_nucl_elec_min,E_nucl_elec)
E_nucl_elec_max = max(E_nucl_elec_max,E_nucl_elec)
END_PROVIDER
BEGIN_PROVIDER [double precision, Eff_pot_mu_elec, (elec_num)]
&BEGIN_PROVIDER [double precision, Eff_pot_mu_elec_simple, (elec_num)]
implicit none
include '../constants.F'
BEGIN_DOC
! Eq.(32)
END_DOC
integer :: i,j
double precision :: rij, mu
mu = mu_erf
Eff_pot_mu_elec = 0.d0
do i=1,elec_num
!DIR$ VECTOR ALIGNED
!DIR$ LOOP COUNT(50)
do j=1,elec_num
rij = elec_dist(j,i)
if(i==j)cycle
Eff_pot_mu_elec(i) = Eff_pot_mu_elec(i) + 0.5d0 * derf(mu * rij) * elec_dist_inv(j,i)
Eff_pot_mu_elec(i) = Eff_pot_mu_elec(i) + 0.5d0 * mu/dsqpi * dexp(-mu*mu*rij*rij)
Eff_pot_mu_elec_simple(i) = Eff_pot_mu_elec(i)
Eff_pot_mu_elec(i) = Eff_pot_mu_elec(i) + 0.5d0 * (- 0.25d0 * (1.d0 - derf(mu*rij))**2.d0 )
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, Eff_pot_mu ]
implicit none
include '../constants.F'
BEGIN_DOC
! Eq.(32)
END_DOC
integer :: i
Eff_pot_mu = 0.d0
do i=1,elec_num
Eff_pot_mu += eff_pot_mu_elec(i)
enddo
Eff_pot_mu_min = min(Eff_pot_mu_min,Eff_pot_mu)
Eff_pot_mu_max = max(Eff_pot_mu_max,Eff_pot_mu)
SOFT_TOUCH Eff_pot_mu_min Eff_pot_mu_max
END_PROVIDER
BEGIN_PROVIDER [double precision, Eff_pot_mu_simple ]
implicit none
BEGIN_DOC
! Eq.(32)
END_DOC
include '../constants.F'
integer :: i
Eff_pot_mu_simple = 0.d0
do i=1,elec_num
Eff_pot_mu_simple += Eff_pot_mu_elec_simple(i)
enddo
Eff_pot_mu_simple_min = min(Eff_pot_mu_simple_min,Eff_pot_mu_simple)
Eff_pot_mu_simple_max = max(Eff_pot_mu_simple_max,Eff_pot_mu_simple)
SOFT_TOUCH Eff_pot_mu_simple_min Eff_pot_mu_simple_max
END_PROVIDER
BEGIN_PROVIDER [double precision, eff_pot_deriv_mu_elec, (elec_num) ]
implicit none
BEGIN_DOC
! Eq.(33)
END_DOC
integer :: i,j
double precision :: rij, mu
mu = mu_erf
eff_pot_deriv_mu_elec = 0.d0
do i = 1, elec_num
do j = 1, elec_num
if(i==j)cycle
rij = elec_dist(i,j)
eff_pot_deriv_mu_elec(i) += 0.5d0 * ( derf(mu * rij) - 1.d0 ) * elec_dist_inv(j,i) &
* ( - elec_dist_vec_x(j,i) * psidet_grad_lapl(1,i) &
- elec_dist_vec_y(j,i) * psidet_grad_lapl(2,i) &
- elec_dist_vec_z(j,i) * psidet_grad_lapl(3,i) ) * psidet_inv
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, three_body_mu ]
implicit none
BEGIN_DOC
! Eq.(30)
END_DOC
integer :: i,j,k
three_body_mu = 0.d0
do i = 1, elec_num
do j = i+1, elec_num
do k = j+1, elec_num
three_body_mu += grad_j_mu_x(i,j) * grad_j_mu_x(i,k)
three_body_mu += grad_j_mu_y(i,j) * grad_j_mu_y(i,k)
three_body_mu += grad_j_mu_z(i,j) * grad_j_mu_z(i,k)
three_body_mu += grad_j_mu_x(j,i) * grad_j_mu_x(j,k)
three_body_mu += grad_j_mu_y(j,i) * grad_j_mu_y(j,k)
three_body_mu += grad_j_mu_z(j,i) * grad_j_mu_z(j,k)
three_body_mu += grad_j_mu_x(k,i) * grad_j_mu_x(k,j)
three_body_mu += grad_j_mu_y(k,i) * grad_j_mu_y(k,j)
three_body_mu += grad_j_mu_z(k,i) * grad_j_mu_z(k,j)
enddo
enddo
enddo
three_body_mu_min = min(three_body_mu_min,three_body_mu)
three_body_mu_max = max(three_body_mu_max,three_body_mu)
SOFT_TOUCH three_body_mu_min three_body_mu_max
END_PROVIDER
BEGIN_PROVIDER [double precision, eff_pot_deriv_mu]
implicit none
BEGIN_DOC
! Eq.(33)
END_DOC
integer :: i
eff_pot_deriv_mu = 0.d0
do i = 1, elec_num
eff_pot_deriv_mu += eff_pot_deriv_mu_elec(i)
enddo
eff_pot_deriv_mu_min = min(eff_pot_deriv_mu_min,eff_pot_deriv_mu)
eff_pot_deriv_mu_max = max(eff_pot_deriv_mu_max,eff_pot_deriv_mu)
SOFT_TOUCH eff_pot_deriv_mu_min eff_pot_deriv_mu_max
END_PROVIDER

1
src/ezfio_interface.irp.f
View File

@ -18,6 +18,7 @@ data = [ \
("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_mu_erf" , "real" , "" ),
("jastrow_jast_a_up_up" , "real" , "" ),
("jastrow_jast_a_up_dn" , "real" , "" ),
("jastrow_jast_b_up_up" , "real" , "" ),

3
src/types.F
View File

@ -11,6 +11,7 @@
integer, parameter :: t_Simple = 21
integer, parameter :: t_None = 22
integer, parameter :: t_Mu = 23
integer, parameter :: t_Core = 24
integer, parameter :: t_Stopped = 0
@ -41,7 +42,7 @@
' ', &
'Simple ', &
'None ', &
' ', &
'Mu ', &
'Core ', &
' '/)