mirror of
https://gitlab.com/scemama/qmcchem.git
synced 2024-09-07 02:24:47 +02:00
315 lines
7.9 KiB
FortranFixed
315 lines
7.9 KiB
FortranFixed
|
BEGIN_PROVIDER [ double precision, elec_mass ]
|
||
|
implicit none
|
||
|
BEGIN_DOC
|
||
|
! Electron mass in the Langevin algorithm
|
||
|
END_DOC
|
||
|
|
||
|
integer :: i
|
||
|
elec_mass=0.d0
|
||
|
do i=1,nucl_num
|
||
|
elec_mass=max(dble(nucl_charge(i)),dble(elec_mass))
|
||
|
enddo
|
||
|
elec_mass=elec_mass**1.5
|
||
|
END_PROVIDER
|
||
|
|
||
|
BEGIN_PROVIDER [ real, elec_impuls_full, (elec_num,4,walk_num)]
|
||
|
implicit none
|
||
|
BEGIN_DOC
|
||
|
! Impulsions used in the Langevin algorithm
|
||
|
END_DOC
|
||
|
integer :: i,k,l
|
||
|
double precision :: temp(elec_num,3)
|
||
|
|
||
|
do k=1,walk_num
|
||
|
call gauss_array(elec_num*3,temp)
|
||
|
do l=1,3
|
||
|
do i=1,elec_num
|
||
|
elec_impuls_full(i,l,k) = time_step*temp(i,l)
|
||
|
enddo
|
||
|
enddo
|
||
|
enddo
|
||
|
|
||
|
END_PROVIDER
|
||
|
|
||
|
BEGIN_PROVIDER [ real, elec_impuls, (elec_num_8,3)]
|
||
|
implicit none
|
||
|
BEGIN_DOC
|
||
|
! Impulsions used in the Langevin algorithm for the current walker
|
||
|
END_DOC
|
||
|
integer :: i,l
|
||
|
do l=1,3
|
||
|
do i=1,elec_num
|
||
|
elec_impuls(i,l) = elec_impuls_full(i,l,walk_i)
|
||
|
enddo
|
||
|
enddo
|
||
|
END_PROVIDER
|
||
|
|
||
|
BEGIN_PROVIDER [ double precision, langevin_sigma, (2) ]
|
||
|
&BEGIN_PROVIDER [ double precision, langevin_c12 ]
|
||
|
&BEGIN_PROVIDER [ double precision, langevin_coef, (4) ]
|
||
|
implicit none
|
||
|
BEGIN_DOC
|
||
|
! Quantities needed for the Langevin algorithm.
|
||
|
END_DOC
|
||
|
langevin_sigma(1) = sqrt (time_step/elec_mass * &
|
||
|
(2.d0- (3.d0-4.d0*time_step_exp + time_step_exp**2)/time_step) )
|
||
|
langevin_sigma(2) = sqrt (elec_mass * (1.d0-time_step_exp**2))
|
||
|
langevin_c12 = (1.d0-time_step_exp)**2/(langevin_sigma(1)*langevin_sigma(2))
|
||
|
langevin_coef(1) = 1.d0/elec_mass*time_step*time_step_exp_sq
|
||
|
langevin_coef(2) = time_step_exp_sq_sq*time_step**2/elec_mass
|
||
|
langevin_coef(3) = langevin_sigma(2)*sqrt(1.d0-langevin_c12**2)
|
||
|
langevin_coef(4) = langevin_c12*langevin_sigma(2)/langevin_sigma(1)
|
||
|
END_PROVIDER
|
||
|
|
||
|
|
||
|
|
||
|
subroutine langevin_step(p,q,accepted,delta_x)
|
||
|
|
||
|
implicit none
|
||
|
|
||
|
double precision, intent(out) :: p,q
|
||
|
logical, intent(out) :: accepted
|
||
|
real, intent(out) :: delta_x
|
||
|
|
||
|
real :: xold(elec_num_8,3)
|
||
|
real :: pold(elec_num_8,3)
|
||
|
double precision :: psiold
|
||
|
double precision :: bold(elec_num_8,3)
|
||
|
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: xold, bold, pold
|
||
|
|
||
|
integer :: i,l
|
||
|
|
||
|
psiold = psi_value
|
||
|
if (psiold == 0.d0) then
|
||
|
call abrt(irp_here,'Walker is on the nodal surface.')
|
||
|
endif
|
||
|
|
||
|
do l=1,3
|
||
|
!DIR$ LOOP COUNT (200)
|
||
|
!DIR$ VECTOR ALIGNED
|
||
|
do i=1,elec_num
|
||
|
xold(i,l) = elec_coord(i,l)
|
||
|
pold(i,l) = elec_impuls(i,l)
|
||
|
enddo
|
||
|
enddo
|
||
|
!DIR$ LOOP COUNT (200)
|
||
|
!DIR$ VECTOR ALIGNED
|
||
|
do i=1,elec_num
|
||
|
bold(i,1) = psi_grad_psi_inv_x(i)
|
||
|
bold(i,2) = psi_grad_psi_inv_y(i)
|
||
|
bold(i,3) = psi_grad_psi_inv_z(i)
|
||
|
enddo
|
||
|
|
||
|
! First move
|
||
|
! W1
|
||
|
call gauss_array(size(xbrown),xbrown)
|
||
|
do l=1,3
|
||
|
!DIR$ VECTOR ALIGNED
|
||
|
!DIR$ LOOP COUNT (200)
|
||
|
do i=1,elec_num
|
||
|
xbrown(i,l) = xbrown(i,l)*langevin_sigma(1)
|
||
|
enddo
|
||
|
enddo
|
||
|
|
||
|
real :: xdiff (elec_num_8,3)
|
||
|
! q(n+1)
|
||
|
delta_x = 0.
|
||
|
do l=1,3
|
||
|
!DIR$ VECTOR ALIGNED
|
||
|
!DIR$ LOOP COUNT (200)
|
||
|
do i=1,elec_num
|
||
|
xdiff(i,l) = pold(i,l)*langevin_coef(1) &
|
||
|
+ bold(i,l)*langevin_coef(2) &
|
||
|
+ xbrown(i,l)
|
||
|
delta_x += xdiff(i,l)*xdiff(i,l)
|
||
|
enddo
|
||
|
enddo
|
||
|
delta_x = sqrt(delta_x)
|
||
|
|
||
|
do l=1,3
|
||
|
!DIR$ VECTOR ALIGNED
|
||
|
!DIR$ LOOP COUNT (200)
|
||
|
do i=1,elec_num
|
||
|
elec_coord(i,l) = elec_coord(i,l) + xdiff(i,l)
|
||
|
enddo
|
||
|
enddo
|
||
|
|
||
|
TOUCH elec_coord
|
||
|
|
||
|
! Second move
|
||
|
! W2
|
||
|
double precision :: gauss
|
||
|
do l=1,3
|
||
|
!DIR$ VECTOR ALIGNED
|
||
|
!DIR$ LOOP COUNT (200)
|
||
|
do i=1,elec_num
|
||
|
xbrown(i,l) = langevin_coef(3)*gauss() + &
|
||
|
langevin_coef(4)*xbrown(i,l)
|
||
|
enddo
|
||
|
enddo
|
||
|
|
||
|
! p(n+1)
|
||
|
!DIR$ LOOP COUNT (200)
|
||
|
do i=1,elec_num
|
||
|
elec_impuls(i,1) = time_step_exp*pold(i,1) &
|
||
|
+ time_step*(bold(i,1)+psi_grad_psi_inv_x(i))*time_step_exp_sq &
|
||
|
+ xbrown(i,1)
|
||
|
elec_impuls(i,2) = time_step_exp*pold(i,2) &
|
||
|
+ time_step*(bold(i,2)+psi_grad_psi_inv_y(i))*time_step_exp_sq &
|
||
|
+ xbrown(i,2)
|
||
|
elec_impuls(i,3) = time_step_exp*pold(i,3) &
|
||
|
+ time_step*(bold(i,3)+psi_grad_psi_inv_z(i))*time_step_exp_sq &
|
||
|
+ xbrown(i,3)
|
||
|
enddo
|
||
|
! TOUCH elec_impuls
|
||
|
! (touch moved after acceptation)
|
||
|
|
||
|
p = 1.d0
|
||
|
|
||
|
double precision :: ratio
|
||
|
ratio = (psi_value/psiold)**2
|
||
|
|
||
|
double precision :: dt_dm
|
||
|
dt_dm = time_step/elec_mass
|
||
|
!DIR$ LOOP COUNT (200)
|
||
|
do i=1,elec_num
|
||
|
|
||
|
double precision :: d11
|
||
|
double precision :: d12
|
||
|
double precision :: temp
|
||
|
double precision :: d21
|
||
|
double precision :: d22
|
||
|
double precision :: re
|
||
|
|
||
|
d11 = -xdiff(i,1) &
|
||
|
+ dt_dm*elec_impuls(i,1)*time_step_exp_sq &
|
||
|
- time_step*dt_dm*psi_grad_psi_inv_x(i)*time_step_exp_sq_sq
|
||
|
|
||
|
d12 = xdiff(i,1) &
|
||
|
- dt_dm*pold(i,1)*time_step_exp_sq &
|
||
|
- time_step*dt_dm*bold(i,1)*time_step_exp_sq_sq
|
||
|
|
||
|
temp = time_step*(psi_grad_psi_inv_x(i)+bold(i,1))*time_step_exp_sq
|
||
|
|
||
|
d21 = -pold(i,1)+elec_impuls(i,1)*time_step_exp - temp
|
||
|
|
||
|
|
||
|
d22 = elec_impuls(i,1)-pold(i,1)*time_step_exp -temp
|
||
|
|
||
|
re = (d11/langevin_sigma(1))**2+(d21/langevin_sigma(2))**2 &
|
||
|
- 2.d0*langevin_c12*d11*d21/(langevin_sigma(1)*langevin_sigma(2))
|
||
|
|
||
|
re -= (d12/langevin_sigma(1))**2+(d22/langevin_sigma(2))**2 &
|
||
|
- 2.d0*langevin_c12*d12*d22/(langevin_sigma(1)*langevin_sigma(2))
|
||
|
|
||
|
re = re / (2.d0*(1.d0-langevin_c12**2))
|
||
|
|
||
|
if (re < 35.d0) then
|
||
|
p=p*exp(-re)
|
||
|
else
|
||
|
p = 0.d0
|
||
|
endif
|
||
|
|
||
|
|
||
|
d11 = -xdiff(i,2) &
|
||
|
+ dt_dm*elec_impuls(i,2)*time_step_exp_sq &
|
||
|
- time_step*dt_dm*psi_grad_psi_inv_y(i)*time_step_exp_sq_sq
|
||
|
|
||
|
d12 = xdiff(i,2) &
|
||
|
- dt_dm*pold(i,2)*time_step_exp_sq &
|
||
|
- time_step*dt_dm*bold(i,2)*time_step_exp_sq_sq
|
||
|
|
||
|
temp = time_step*(psi_grad_psi_inv_y(i)+bold(i,2))*time_step_exp_sq
|
||
|
|
||
|
d21 = -pold(i,2)+elec_impuls(i,2)*time_step_exp - temp
|
||
|
|
||
|
|
||
|
d22 = elec_impuls(i,2)-pold(i,2)*time_step_exp -temp
|
||
|
|
||
|
re = (d11/langevin_sigma(1))**2+(d21/langevin_sigma(2))**2 &
|
||
|
- 2.d0*langevin_c12*d11*d21/(langevin_sigma(1)*langevin_sigma(2))
|
||
|
|
||
|
re -= (d12/langevin_sigma(1))**2+(d22/langevin_sigma(2))**2 &
|
||
|
- 2.d0*langevin_c12*d12*d22/(langevin_sigma(1)*langevin_sigma(2))
|
||
|
|
||
|
re = re / (2.d0*(1.d0-langevin_c12**2))
|
||
|
|
||
|
if (re < 35.d0) then
|
||
|
p=p*exp(-re)
|
||
|
else
|
||
|
p = 0.d0
|
||
|
endif
|
||
|
|
||
|
d11 = -xdiff(i,3) &
|
||
|
+ dt_dm*elec_impuls(i,3)*time_step_exp_sq &
|
||
|
- time_step*dt_dm*psi_grad_psi_inv_z(i)*time_step_exp_sq_sq
|
||
|
|
||
|
d12 = xdiff(i,3) &
|
||
|
- dt_dm*pold(i,3)*time_step_exp_sq &
|
||
|
- time_step*dt_dm*bold(i,3)*time_step_exp_sq_sq
|
||
|
|
||
|
temp = time_step*(psi_grad_psi_inv_z(i)+bold(i,3))*time_step_exp_sq
|
||
|
|
||
|
d21 = -pold(i,3)+elec_impuls(i,3)*time_step_exp - temp
|
||
|
|
||
|
|
||
|
d22 = elec_impuls(i,3)-pold(i,3)*time_step_exp -temp
|
||
|
|
||
|
re = (d11/langevin_sigma(1))**2+(d21/langevin_sigma(2))**2 &
|
||
|
- 2.d0*langevin_c12*d11*d21/(langevin_sigma(1)*langevin_sigma(2))
|
||
|
|
||
|
re -= (d12/langevin_sigma(1))**2+(d22/langevin_sigma(2))**2 &
|
||
|
- 2.d0*langevin_c12*d12*d22/(langevin_sigma(1)*langevin_sigma(2))
|
||
|
|
||
|
re = re / (2.d0*(1.d0-langevin_c12**2))
|
||
|
|
||
|
if (re < 35.d0) then
|
||
|
p=p*exp(-re)
|
||
|
else
|
||
|
p = 0.d0
|
||
|
endif
|
||
|
|
||
|
|
||
|
double precision :: sumup, sumdn
|
||
|
sumup = elec_impuls(i,1)*elec_impuls(i,1) &
|
||
|
+ elec_impuls(i,2)*elec_impuls(i,2) &
|
||
|
+ elec_impuls(i,3)*elec_impuls(i,3)
|
||
|
sumdn = pold(i,1)*pold(i,1) + pold(i,2)*pold(i,2) + pold(i,3)*pold(i,3)
|
||
|
re = exp(-(sumup-sumdn)/(2.d0*elec_mass))
|
||
|
p = p*re
|
||
|
|
||
|
enddo
|
||
|
|
||
|
|
||
|
p = min (1.d0,ratio*p)
|
||
|
q = 1.d0 - p
|
||
|
|
||
|
double precision :: qmc_ranf
|
||
|
accepted = p > qmc_ranf()
|
||
|
if (accepted) then
|
||
|
accepted_num += 1.
|
||
|
SOFT_TOUCH elec_impuls
|
||
|
else
|
||
|
do l=1,3
|
||
|
!DIR$ VECTOR ALIGNED
|
||
|
!DIR$ LOOP COUNT(200)
|
||
|
do i=1,elec_num
|
||
|
elec_coord(i,l) = xold(i,l)
|
||
|
elec_impuls(i,l) = -pold(i,l)
|
||
|
enddo
|
||
|
enddo
|
||
|
!DIR$ VECTOR ALIGNED
|
||
|
!DIR$ LOOP COUNT(200)
|
||
|
do i=1,elec_num
|
||
|
psi_grad_psi_inv_x(i) = bold(i,1)
|
||
|
psi_grad_psi_inv_y(i) = bold(i,2)
|
||
|
psi_grad_psi_inv_z(i) = bold(i,3)
|
||
|
enddo
|
||
|
psi_value = psiold
|
||
|
rejected_num += 1.
|
||
|
SOFT_TOUCH elec_coord psi_grad_psi_inv_x psi_grad_psi_inv_y psi_grad_psi_inv_z psi_value
|
||
|
endif
|
||
|
|
||
|
end
|
||
|
|