Dice Fortran Backend Documentation
quantum_solvers.f90
Go to the documentation of this file.
1 ! ---------------------------------------------------------------------------------------
2 ! Dice Quantum Monte Carlo
3 ! ---------------------------------------------------------------------------------------
4 ! MODULE: solvers
5 !
6 ! DESCRIPTION:
9 !
14 ! ---------------------------------------------------------------------------------------
15 
16 module solvers
17 
18  use constants
19  use shared_data
20 
21  implicit none
22 
23 contains
24 
25 ! ---------------------------------------------------------------------------------------
26 ! SECTION: QHO Solvers
27 ! ---------------------------------------------------------------------------------------
28 
29 
30  ! -----------------------------------------------------------------------------------
31  ! ROUTINE: QHO_Prob
32  !
33  ! DESCRIPTION:
36  !
40  !
41  ! PARAMETERS:
46  ! -----------------------------------------------------------------------------------
47  function qho_prob(alpha, x_current, x_next) result(Prob)
48 
49  implicit none
50 
51  real(real64), intent(in) :: alpha
52  real(real64), intent(in) :: x_current, x_next
53 
54  real(real64) :: prob
55 
56  prob = exp(-2*alpha*(x_next*x_next - x_current*x_current))
57 
58  end function qho_prob
59 
60  ! -----------------------------------------------------------------------------------
61  ! ROUTINE: QHO_Energy
62  !
63  ! DESCRIPTION:
66  !
70  !
71  ! PARAMETERS:
75  ! -----------------------------------------------------------------------------------
76  function qho_energy(alpha, x) result(E_loc)
77 
78  implicit none
79 
80  real(real64), intent(in) :: alpha
81  real(real64), intent(in) :: x
82 
83  real(real64) :: e_loc
84 
85  e_loc = alpha + (x*x)*(0.5 - (2*alpha*alpha))
86 
87  end function qho_energy
88 
89 ! ---------------------------------------------------------------------------------------
90 ! SECTION: Hydrogenic Wavefunctions
91 ! ---------------------------------------------------------------------------------------
92 
93  ! -----------------------------------------------------------------------------------
94  ! ROUTINE: H1s_wfn
95  !
96  ! DESCRIPTION:
99  !
105  !
106  ! PARAMETERS:
110  ! -----------------------------------------------------------------------------------
111  function h1s_wfn(r_elec, R_nuc)
112 
113  implicit none
114 
115  real(real64), dimension(3), intent(in) :: r_elec, r_nuc
116 
117  real(real64) :: h1s_wfn, expfac
118 
119  expfac = -norm2(r_elec - r_nuc)
120 
121  h1s_wfn = pi_dash * exp(expfac)
122 
123  end function h1s_wfn
124 
125  ! -----------------------------------------------------------------------------------
126  ! ROUTINE: H2s_wfn
127  !
128  ! DESCRIPTION:
131  !
143  !
144  ! PARAMETERS:
148  ! -----------------------------------------------------------------------------------
149  function h2s_wfn(r_elec, R_nuc)
150 
151  implicit none
152 
153  real(real64), dimension(3), intent(in) :: r_elec, r_nuc
154 
155  real(real64) :: h2s_wfn, prefac, expfac, r_norm
156 
157  r_norm = norm2(r_elec - r_nuc)
158  prefac = pi_dash_2 * (2 - r_norm)
159  expfac = -(r_norm / 2)
160 
161  h2s_wfn = prefac * exp(expfac)
162 
163  end function h2s_wfn
164 
165 ! ---------------------------------------------------------------------------------------
166 ! SECTION: H_2_plus Solvers
167 ! ---------------------------------------------------------------------------------------
168 
169  ! -----------------------------------------------------------------------------------
170  ! ROUTINE: H_2_plus_wfn
171  !
172  ! DESCRIPTION:
175  !
182  !
183  ! PARAMETERS:
189  ! -----------------------------------------------------------------------------------
190  function h_2_plus_wfn(c, r, R_a, R_b) result(wfn)
191 
192  implicit none
193 
194  real(real64), intent(in) :: c
195  real(real64), dimension(3), intent(in) :: r
196  real(real64), dimension(3), intent(in) :: r_a, r_b
197 
198  real(real64) :: wfn
199 
200  wfn = (c * h1s_wfn(r, r_a)) + (sqrt(1 - (c*c)) * h1s_wfn(r, r_b))
201 
202  end function h_2_plus_wfn
203 
204  ! -----------------------------------------------------------------------------------
205  ! ROUTINE: H_2_plus_Prob
206  !
207  ! DESCRIPTION:
210  !
216  !
217  ! PARAMETERS:
224  ! -----------------------------------------------------------------------------------
225  function h_2_plus_prob(c, r_current, r_next, R_a, R_b) result(Prob)
226 
227  implicit none
228 
229  real(real64), intent(in) :: c
230  real(real64), dimension(3), intent(in) :: r_current, r_next
231  real(real64), dimension(3), intent(in) :: r_a, r_b
232 
233  real(real64) :: prob
234 
235  prob = (h_2_plus_wfn(c, r_next, r_a, r_b) / &
236  h_2_plus_wfn(c, r_current, r_a, r_b))**2
237 
238  end function h_2_plus_prob
239 
240  ! -----------------------------------------------------------------------------------
241  ! ROUTINE: H_2_plus_Energy
242  !
243  ! DESCRIPTION:
246  !
268  !
269  ! PARAMETERS:
276  ! -----------------------------------------------------------------------------------
277  subroutine h_2_plus_energy(c, r, R_a, R_b, E_loc, gradient)
278 
279  implicit none
280 
281  real(real64), intent(in) :: c
282  real(real64), dimension(3), intent(in) :: r
283  real(real64), dimension(3), intent(in) :: R_a, R_b
284 
285  real(real64), intent(out) :: E_loc
286  real(real64), intent(out) :: gradient
287 
288  real(real64), dimension(2, 3) :: r_diffs ! Wavefunction finite differences wrt position
289  real(real64), dimension(2) :: c_diffs ! Wavefunction finite differences wrt parameter
290  real(real64), dimension(3) :: r_temp ! Temporary position for finite differences
291  real(real64) :: laplacian ! 2nd derivative of trial wavefunction
292  real(real64) :: coulomb ! Coulomb terms of Hamiltonian
293  real(real64) :: c_temp ! Temporary parameter value for finite differences
294  real(real64) :: wfn, wfnx2 ! Wavefunction at current walker position
295  integer :: i ! Loop integer
296 
297  ! Calculate trial wavefunction at provided point.
298  wfn = h_2_plus_wfn(c, r, r_a, r_b)
299  wfnx2 = 2 * wfn
300 
301  ! Calculate trial wavefunction in finite difference directions.
302  do i = 1, 3
303  r_temp(:) = r(:)
304  r_temp(i) = r(i) - h2plus%ham_fdstep
305  r_diffs(1, i) = h_2_plus_wfn(c, r_temp, r_a, r_b)
306 
307  r_temp(i) = r(i) + h2plus%ham_fdstep
308  r_diffs(2, i) = h_2_plus_wfn(c, r_temp, r_a, r_b)
309  end do
310 
311  ! Assemble Hamiltonian
312  laplacian = ((r_diffs(1, 1) - wfnx2 + r_diffs(2, 1)) + &
313  (r_diffs(1, 2) - wfnx2 + r_diffs(2, 2)) + &
314  (r_diffs(1, 3) - wfnx2 + r_diffs(2, 3))) / (h2plus%ham_fdstep**2)
315  laplacian = laplacian * (-0.5_real64)
316 
317  coulomb = -(1 / norm2(r - r_a))*wfn - (1 / norm2(r - r_b))*wfn
318 
319  e_loc = ((laplacian + coulomb) / wfn) + (1 / norm2(r_a - r_b))
320 
321  ! Check if gradient is needed
322  if (h2plus%auto_params) then
323  ! Calculate wavefunction gradient wrt c by central difference.
324  c_temp = c - h2plus%c_fdstep
325  c_diffs(1) = h_2_plus_wfn(c_temp, r, r_a, r_b)
326 
327  c_temp = c + h2plus%c_fdstep
328  c_diffs(2) = h_2_plus_wfn(c_temp, r, r_a, r_b)
329 
330  gradient = (c_diffs(2) - c_diffs(1)) / (2.0_real64 * h2plus%c_fdstep)
331  else
332  gradient = 0.0_real64
333  end if
334 
335  end subroutine h_2_plus_energy
336 
337  ! -----------------------------------------------------------------------------------
338  ! ROUTINE: H_2_plus_update_c
339  !
340  ! DESCRIPTION:
344  !
348  !
349  ! PARAMETERS:
353  ! -----------------------------------------------------------------------------------
354  subroutine h_2_plus_update_c(c_old, gradient, c_new)
355 
356  implicit none
357 
358  real(real64), intent(in) :: c_old
359  real(real64), intent(in) :: gradient
360  real(real64), intent(out) :: c_new
361 
362  real(real64) :: damping = 1.0_real64
363 
364  ! Update c by dampened steepest descent.
365  c_new = c_old - (damping * gradient)
366 
367  end subroutine h_2_plus_update_c
368 
369 ! ---------------------------------------------------------------------------------------
370 ! SECTION: H_2 Solvers
371 ! ---------------------------------------------------------------------------------------
372 
373  ! -----------------------------------------------------------------------------------
374  ! ROUTINE: H_2_cusp
375  !
376  ! DESCRIPTION:
380  !
389  !
390  ! PARAMETERS:
393  ! -----------------------------------------------------------------------------------
394  function h_2_cusp(s) result(a)
395 
396  implicit none
397 
398  real(real64), intent(in) :: s
399 
400  real(real64) :: a
401 
402  real(real64) :: f ! Value of function at current iteration
403  real(real64) :: f_dash ! Value of derivative at current iteration
404  real(real64) :: emsoa ! e^(-s/a)
405  real(real64) :: tol ! Tolerance for stopping iteration
406  integer :: maxiter ! Maximum number of iterations
407  integer :: k ! Loop integer
408 
409  a = 0.5_real64 ! guess based on s=0 solution
410  maxiter = 100
411  tol = 1e-5_real64
412 
413  if (i_log) then
414  write(logid, "('Starting Newton-Raphson search for parameter a...')")
415  write(logid, *) 'a abs(f)'
416  end if
417 
418  ! Calculate initial value of function based on guess
419  f = (1.0_real64 / (1.0_real64 + exp(-s/a))) - a
420  ! Iterate until converged
421  do k = 1, maxiter
422 
423  if (abs(f) .lt. tol) then
424  if (i_log) then
425  write(logid, "('Value of parameter a converged.')")
426  write(logid, "('')")
427  end if
428  exit
429  end if
430 
431  ! Calculate new values
432  emsoa = exp(-s/a)
433  f = (1.0_real64 / (1.0_real64 + emsoa)) - a
434  f_dash = (-s * emsoa) / ((1.0_real64 + emsoa*emsoa * a*a)) - 1
435 
436  ! Update value of a
437  a = a - (f / f_dash)
438 
439  if (i_log) then
440  write(logid, "(F12.8, ' ', F12.8)") a, abs(f)
441  end if
442 
443  end do
444 
445  ! Crash out if not converged
446  if (k .ge. maxiter) then
447  error stop "Newton-Raphson search for parameter a has failed to converge."
448  end if
449 
450  end function h_2_cusp
451 
452  ! -----------------------------------------------------------------------------------
453  ! ROUTINE: H_2_wfn
454  !
455  ! DESCRIPTION:
458  !
478  !
479  ! PARAMETERS:
487  ! -----------------------------------------------------------------------------------
488  function h_2_wfn(a, beta, r_1, r_2, R_a, R_b) result(wfn)
489 
490  implicit none
491 
492  real(real64), intent(in) :: a
493  real(real64), intent(in) :: beta
494  real(real64), dimension(3), intent(in) :: r_1, r_2
495  real(real64), dimension(3), intent(in) :: r_a, r_b
496 
497  real(real64) :: wfn
498 
499  real(real64) :: s ! Bond length
500  real(real64) :: alpha ! Wavefunction parameter, set to 2
501  real(real64) :: phi_r1, phi_r2, jastrow ! Wavefunction components
502 
503  s = norm2(r_a - r_b)
504  alpha = 2.0_real64
505 
506  phi_r1 = exp(-norm2(r_1 - r_a) / a) + exp(-norm2(r_1 - r_b) / a)
507  phi_r2 = exp(-norm2(r_2 - r_a) / a) + exp(-norm2(r_2 - r_b) / a)
508 
509  jastrow = exp(norm2(r_1 - r_2) / (alpha * (1.0_real64 + beta*norm2(r_1 - r_2))))
510 
511  wfn = phi_r1 * phi_r2 * jastrow
512 
513  end function h_2_wfn
514 
515  ! -----------------------------------------------------------------------------------
516  ! ROUTINE: H_2_Prob
517  !
518  ! DESCRIPTION:
521  !
527  !
528  ! PARAMETERS:
538  ! -----------------------------------------------------------------------------------
539  function h_2_prob(a, beta, r_1_current, r_2_current, r_1_next, r_2_next, &
540  R_a, R_b) result(Prob)
541 
542  implicit none
543 
544  real(real64), intent(in) :: a
545  real(real64), intent(in) :: beta
546  real(real64), dimension(3), intent(in) :: r_1_current, r_2_current
547  real(real64), dimension(3), intent(in) :: r_1_next, r_2_next
548  real(real64), dimension(3), intent(in) :: r_a, r_b
549 
550  real(real64) :: prob
551 
552  prob = (h_2_wfn(a, beta, r_1_next, r_2_next, r_a, r_b) / &
553  h_2_wfn(a, beta, r_1_current, r_2_current, r_a, r_b))**2
554 
555  end function h_2_prob
556 
557  ! -----------------------------------------------------------------------------------
558  ! ROUTINE: H_2_Energy
559  !
560  ! DESCRIPTION:
563  !
580  !
581  ! PARAMETERS:
590  ! -----------------------------------------------------------------------------------
591  subroutine h_2_energy(a, beta, r_1, r_2, R_a, R_b, E_loc, DphiDbeta)
592 
593  implicit none
594 
595  real(real64), intent(in) :: a
596  real(real64), intent(in) :: beta
597  real(real64), dimension(3), intent(in) :: r_1, r_2
598  real(real64), dimension(3), intent(in) :: R_a, R_b
599 
600  real(real64), intent(out) :: E_loc
601  real(real64), intent(out) :: DphiDbeta
602 
603  real(real64) :: E1, E2, E3, E4, E5, E6, E7, E8
604  real(real64) :: phi_1a, phi_1b, phi_1
605  real(real64) :: phi_2a, phi_2b, phi_2
606  real(real64) :: r_1a, r_1b, r_2a, r_2b, r_12
607  real(real64), dimension(3) :: r_1a_hat, r_1b_hat, r_2a_hat, r_2b_hat, r_12_hat
608 
609  r_1a = norm2(r_1 - r_a)
610  r_1b = norm2(r_1 - r_b)
611  r_2a = norm2(r_2 - r_a)
612  r_2b = norm2(r_2 - r_b)
613  r_12 = norm2(r_1 - r_2)
614 
615  r_1a_hat = (r_1 - r_a) / spread(r_1a, 1, 3)
616  r_1b_hat = (r_1 - r_b) / spread(r_1b, 1, 3)
617  r_2a_hat = (r_2 - r_a) / spread(r_2a, 1, 3)
618  r_2b_hat = (r_2 - r_b) / spread(r_2b, 1, 3)
619  r_12_hat = (r_1 - r_2) / spread(r_12, 1, 3)
620 
621  phi_1a = exp(-r_1a / a)
622  phi_1b = exp(-r_1b / a)
623  phi_2a = exp(-r_2a / a)
624  phi_2b = exp(-r_2b / a)
625 
626  phi_1 = phi_1a + phi_1b
627  phi_2 = phi_2a + phi_2b
628 
629  e1 = 1.0_real64 / (a*a)
630  e2 = ((phi_1a / r_1a) + (phi_1b / r_1b)) / (a * phi_1)
631  e3 = ((phi_2a / r_2a) + (phi_2b / r_2b)) / (a * phi_2)
632  e4 = (1.0_real64 / r_1a) + (1.0_real64 / r_1b) + (1.0_real64 / r_2a) + (1.0_real64 / r_2b)
633  e5 = 1.0_real64 / r_12
634  e6 = ((((phi_1a * sum(r_1a_hat * r_12_hat, dim=1)) + (phi_1b * sum(r_1b_hat * r_12_hat, dim=1))) / phi_1) &
635  - (((phi_2a * sum(r_2a_hat * r_12_hat, dim=1)) + (phi_2b * sum(r_2b_hat * r_12_hat, dim=1))) / phi_2)) &
636  / (2.0_real64 * a * (1.0_real64 + beta*r_12)**2)
637  e7 = (((4.0_real64*beta + 1.0_real64) * r_12) + 4.0_real64) &
638  / (4.0_real64 * r_12 * (1.0_real64 + beta*r_12)**4)
639  e8 = (1.0_real64 / norm2(r_a - r_b))
640 
641  e_loc = -e1 + e2 + e3 - e4 + e5 + e6 - e7 + e8
642 
643  dphidbeta = -(r_12*r_12) / (2.0_real64*(1.0_real64 + beta*r_12)**2)
644 
645  end subroutine h_2_energy
646 
647  ! -----------------------------------------------------------------------------------
648  ! ROUTINE: H_2_update_beta
649  !
650  ! DESCRIPTION:
654  !
658  !
659  ! PARAMETERS:
663  ! -----------------------------------------------------------------------------------
664  subroutine h_2_update_beta(beta_old, gradient, beta_new)
665 
666  implicit none
667 
668  real(real64), intent(in) :: beta_old
669  real(real64), intent(in) :: gradient
670  real(real64), intent(out) :: beta_new
671 
672  real(real64) :: damping = 1.0_real64
673 
674  beta_new = beta_old - (damping * gradient)
675 
676  end subroutine h_2_update_beta
677 
678 end module solvers
constants
Contains mathematical constants needed throughout the code, all of which are globally accessible.
Definition: shared_data.f90:16
constants::pi_dash_2
real(real64), parameter pi_dash_2
1 / (4*sqrt(2*pi)) - multiplication factor used by H2s_wfn
Definition: shared_data.f90:26
constants::pi_dash
real(real64), parameter pi_dash
1 / sqrt(pi) - multiplication factor used by H1s_wfn
Definition: shared_data.f90:25
shared_data
Contains derived types, and global variables to store input values of simulation parameters.
Definition: shared_data.f90:49
shared_data::h2plus
type(h2plus_type) h2plus
To store inputs for H2 ion system.
Definition: shared_data.f90:139
shared_data::i_log
logical i_log
Indicates whether the processor needs to write the logfile.
Definition: shared_data.f90:161
shared_data::logid
integer logid
To store logfile unit.
Definition: shared_data.f90:160
solvers
Contains local energy and transition probability solvers for available problems.
Definition: quantum_solvers.f90:16
solvers::h_2_cusp
real(real64) function h_2_cusp(s)
Calculates a value for a given the bond length s by Newton-Raphson iteration.
Definition: quantum_solvers.f90:395
solvers::h_2_update_beta
subroutine h_2_update_beta(beta_old, gradient, beta_new)
Updates the parameter in the problem by dampened steepest descent.
Definition: quantum_solvers.f90:665
solvers::qho_energy
real(real64) function qho_energy(alpha, x)
Calculates local energy for quantum harmonic oscillator problem.
Definition: quantum_solvers.f90:77
solvers::h_2_plus_prob
real(real64) function h_2_plus_prob(c, r_current, r_next, R_a, R_b)
Calculates transition probability for the problem.
Definition: quantum_solvers.f90:226
solvers::h_2_plus_wfn
real(real64) function h_2_plus_wfn(c, r, R_a, R_b)
Calculates trial wavefunction for the problem.
Definition: quantum_solvers.f90:191
solvers::qho_prob
real(real64) function qho_prob(alpha, x_current, x_next)
Calculates transition probability for quantum harmonic oscillator problem.
Definition: quantum_solvers.f90:48
solvers::h1s_wfn
real(real64) function h1s_wfn(r_elec, R_nuc)
Calculates value of the wavefunction of a hydrogen 1s orbital.
Definition: quantum_solvers.f90:112
solvers::h2s_wfn
real(real64) function h2s_wfn(r_elec, R_nuc)
[DEPRECATED] Calculates value of the wavefunction of a hydrogen 2s orbital.
Definition: quantum_solvers.f90:150
solvers::h_2_energy
subroutine h_2_energy(a, beta, r_1, r_2, R_a, R_b, E_loc, DphiDbeta)
{UPDATE ME} Calculates local energy for the problem.
Definition: quantum_solvers.f90:592
solvers::h_2_plus_energy
subroutine h_2_plus_energy(c, r, R_a, R_b, E_loc, gradient)
Calculates local energy for the problem.
Definition: quantum_solvers.f90:278
solvers::h_2_prob
real(real64) function h_2_prob(a, beta, r_1_current, r_2_current, r_1_next, r_2_next, R_a, R_b)
Calculates transition probability for the problem.
Definition: quantum_solvers.f90:541
solvers::h_2_wfn
real(real64) function h_2_wfn(a, beta, r_1, r_2, R_a, R_b)
Calculates trial wavefunction for the problem.
Definition: quantum_solvers.f90:489
solvers::h_2_plus_update_c
subroutine h_2_plus_update_c(c_old, gradient, c_new)
Updates the parameter in the problem by dampened steepest descent.
Definition: quantum_solvers.f90:355