LCPQ/qmcchem
10
1
Fork 1
mirror of https://gitlab.com/scemama/qmcchem.git synced 2025-03-25 10:06:32 +01:00
Code Releases Activity
qmcchem/src/SAMPLING/langevin_step.irp.f

315 lines
7.9 KiB
Fortran
Raw Permalink Blame History

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
View Git Blame Copy Permalink