mirror of
https://github.com/TREX-CoE/qmckl.git
synced 2025-04-29 20:04:50 +02:00
optimizations
This commit is contained in:
Emiel Slootman
parent
cf8eb9ffea
commit
040b7a3fe9
@ -8087,6 +8087,8 @@ qmckl_exit_code qmckl_provide_ao_basis_ao_hessian(qmckl_context context)
|
||||
return qmckl_failwith( context, rc, "qmckl_provide_ao_basis_shell_hessian", NULL);
|
||||
}
|
||||
|
||||
memset(ctx->ao_basis.ao_hessian, 0, mem_info.size);
|
||||
|
||||
rc = qmckl_compute_ao_hessian_doc(context,
|
||||
ctx->ao_basis.ao_num,
|
||||
ctx->ao_basis.shell_num,
|
||||
@ -8206,11 +8208,13 @@ function qmckl_compute_ao_hessian_doc(context, &
|
||||
! Don't compute polynomials when the radial part is zero.
|
||||
<<fortran_cutoff>>
|
||||
|
||||
!ao_hessian = 0.d0
|
||||
|
||||
do ipoint = 1, point_num
|
||||
e_coord(1) = coord(ipoint,1)
|
||||
e_coord(2) = coord(ipoint,2)
|
||||
e_coord(3) = coord(ipoint,3)
|
||||
ao_hessian(:,:,ipoint,:,:) = 0.d0
|
||||
!o_hessian(:,:,ipoint,:,:) = 0.d0
|
||||
do inucl=1,nucl_num
|
||||
n_coord(1) = nucl_coord(inucl,1)
|
||||
n_coord(2) = nucl_coord(inucl,2)
|
||||
@ -8223,6 +8227,7 @@ function qmckl_compute_ao_hessian_doc(context, &
|
||||
|
||||
r2 = x*x + y*y + z*z
|
||||
if (r2 > cutoff*nucleus_range(inucl)) then
|
||||
!ao_hessian(:,:,ipoint,:,inucl) = 0.d0
|
||||
cycle
|
||||
end if
|
||||
|
||||
@ -8266,24 +8271,6 @@ function qmckl_compute_ao_hessian_doc(context, &
|
||||
poly_vgl(4,il) * shell_hessian(ishell,3,i,ipoint) &
|
||||
)) * ao_factor(k)
|
||||
|
||||
! only p gggp wrong
|
||||
!ao_hessian(k,4,ipoint,i,inucl) = (&
|
||||
! poly_vgl(1,il) * shell_hessian(ishell,4,i,ipoint) ) * ao_factor(k)
|
||||
|
||||
!ao_hessian(k,4,ipoint,i,inucl) = (&
|
||||
! poly_vgl(5,il) * shell_vgl(ishell,i+1,ipoint) + &
|
||||
! poly_vgl(i+1,il) * shell_vgl(ishell,5,ipoint)) * ao_factor(k)
|
||||
|
||||
!ao_hessian(k,4,ipoint,i,inucl) = (&
|
||||
! 2.d0 * (&
|
||||
! poly_hessian(1,i,il) * shell_vgl(ishell,2,ipoint) + &
|
||||
! poly_vgl(2,il) * shell_hessian(ishell,1,i,ipoint) + &
|
||||
! poly_hessian(2,i,il) * shell_vgl(ishell,3,ipoint) + &
|
||||
! poly_vgl(3,il) * shell_hessian(ishell,2,i,ipoint) + &
|
||||
! poly_hessian(3,i,il) * shell_vgl(ishell,4,ipoint) + &
|
||||
! poly_vgl(4,il) * shell_hessian(ishell,3,i,ipoint) &
|
||||
! )) * ao_factor(k)
|
||||
|
||||
end do
|
||||
|
||||
k = k+1
|
||||
|
||||
@ -35,6 +35,7 @@ int main() {
|
||||
#include "qmckl_mo_private_type.h"
|
||||
#include "qmckl_jastrow_champ_private_type.h"
|
||||
#include "qmckl_jastrow_champ_single_private_type.h"
|
||||
#include "qmckl_jastrow_champ_quad_private_type.h"
|
||||
#include "qmckl_forces_private_type.h"
|
||||
#include "qmckl_determinant_private_type.h"
|
||||
#include "qmckl_local_energy_private_type.h"
|
||||
@ -45,6 +46,7 @@ int main() {
|
||||
#include "qmckl_mo_private_func.h"
|
||||
#include "qmckl_jastrow_champ_private_func.h"
|
||||
#include "qmckl_jastrow_champ_single_private_func.h"
|
||||
#include "qmckl_jastrow_champ_quad_private_func.h"
|
||||
#include "qmckl_forces_private_func.h"
|
||||
#include "qmckl_determinant_private_func.h"
|
||||
#include "qmckl_local_energy_private_func.h"
|
||||
@ -134,6 +136,7 @@ typedef struct qmckl_context_struct {
|
||||
/* Points */
|
||||
qmckl_point_struct point;
|
||||
qmckl_jastrow_champ_single_struct single_point;
|
||||
qmckl_jastrow_champ_quad_struct quad_point;
|
||||
|
||||
/* -- Molecular system -- */
|
||||
qmckl_nucleus_struct nucleus;
|
||||
|
||||
@ -4342,7 +4342,7 @@ integer(qmckl_exit_code) function qmckl_compute_forces_jastrow_delta_p_doc( &
|
||||
|
||||
|
||||
integer*8 :: i, a, j, l, k, p, m, n, nw, num
|
||||
double precision :: tmp
|
||||
double precision :: tmp, accu
|
||||
integer*8 :: LDA, LDB, LDC
|
||||
|
||||
num = num_in + 1
|
||||
@ -4378,9 +4378,11 @@ integer(qmckl_exit_code) function qmckl_compute_forces_jastrow_delta_p_doc( &
|
||||
forces_delta_p(j,k,a,l,m,nw) = delta_e(j) * (-een_rescaled_single_n_gl(k,a,l,nw)) - &
|
||||
een_rescaled_e(j,num,m,nw) * delta_n_gl
|
||||
end do
|
||||
accu = 0.d0
|
||||
do j = 1, elec_num
|
||||
forces_delta_p(num,k,a,l,m,nw) = forces_delta_p(num,k,a,l,m,nw) - delta_e(j) * een_rescaled_n_gl(j,k,a,l,nw)
|
||||
accu = accu - delta_e(j) * een_rescaled_n_gl(j,k,a,l,nw)
|
||||
end do
|
||||
forces_delta_p(num,k,a,l,m,nw) = forces_delta_p(num,k,a,l,m,nw) + accu
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
@ -5342,6 +5344,8 @@ qmckl_exit_code qmckl_provide_forces_ao_value(qmckl_context context)
|
||||
return qmckl_failwith( context, rc, "qmckl_provide_ao_basis_ao_vgl", NULL);
|
||||
}
|
||||
|
||||
memset(ctx->forces.forces_ao_value, 0, mem_info.size);
|
||||
|
||||
rc = qmckl_compute_forces_ao_value_doc(context,
|
||||
ctx->ao_basis.ao_num,
|
||||
ctx->ao_basis.shell_num,
|
||||
@ -5419,7 +5423,7 @@ function qmckl_compute_forces_ao_value_doc(context, &
|
||||
integer , allocatable :: ao_index(:)
|
||||
allocate(ao_index(ao_num))
|
||||
|
||||
forces_ao_value = 0.d0
|
||||
!forces_ao_value = 0.d0
|
||||
|
||||
! Pre-computed data
|
||||
do l=0,20
|
||||
@ -5636,6 +5640,71 @@ qmckl_get_forces_mo_value(qmckl_context context,
|
||||
end interface
|
||||
#+end_src
|
||||
|
||||
|
||||
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
|
||||
qmckl_exit_code
|
||||
qmckl_get_forces_mo_value_inplace (qmckl_context context,
|
||||
double* const forces_mo_value,
|
||||
const int64_t size_max);
|
||||
#+end_src
|
||||
|
||||
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
||||
qmckl_exit_code
|
||||
qmckl_get_forces_mo_value_inplace (qmckl_context context,
|
||||
double* const forces_mo_value,
|
||||
const int64_t size_max)
|
||||
{
|
||||
|
||||
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
||||
return qmckl_failwith( context,
|
||||
QMCKL_INVALID_CONTEXT,
|
||||
"qmckl_get_forces_mo_value",
|
||||
NULL);
|
||||
}
|
||||
|
||||
qmckl_exit_code rc;
|
||||
|
||||
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
|
||||
assert (ctx != NULL);
|
||||
|
||||
const int64_t sze = ctx->point.num * 3 * ctx->mo_basis.mo_num * ctx->nucleus.num;
|
||||
if (size_max < sze) {
|
||||
return qmckl_failwith( context,
|
||||
QMCKL_INVALID_ARG_3,
|
||||
"qmckl_get_mo_basis_mo_value",
|
||||
"input array too small");
|
||||
}
|
||||
|
||||
rc = qmckl_context_touch(context);
|
||||
if (rc != QMCKL_SUCCESS) return rc;
|
||||
|
||||
double* old_array = ctx->forces.forces_mo_value;
|
||||
|
||||
ctx->forces.forces_mo_value = forces_mo_value;
|
||||
|
||||
rc = qmckl_provide_forces_mo_value(context);
|
||||
if (rc != QMCKL_SUCCESS) return rc;
|
||||
|
||||
ctx->forces.forces_mo_value = old_array;
|
||||
|
||||
return QMCKL_SUCCESS;
|
||||
}
|
||||
#+end_src
|
||||
|
||||
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
|
||||
interface
|
||||
integer(qmckl_exit_code) function qmckl_get_forces_mo_value_inplace (context, &
|
||||
forces_mo_value, size_max) bind(C)
|
||||
use, intrinsic :: iso_c_binding
|
||||
import
|
||||
implicit none
|
||||
integer (c_int64_t) , intent(in) , value :: context
|
||||
real(c_double), intent(out) :: forces_mo_value(*)
|
||||
integer (c_int64_t) , intent(in) , value :: size_max
|
||||
end function qmckl_get_forces_mo_value_inplace
|
||||
end interface
|
||||
#+end_src
|
||||
|
||||
** Provide
|
||||
|
||||
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
|
||||
@ -5683,6 +5752,7 @@ qmckl_exit_code qmckl_provide_forces_mo_value(qmckl_context context)
|
||||
if (ctx->forces.forces_mo_value == NULL) {
|
||||
|
||||
double* forces_mo_value = (double*) qmckl_malloc(context, mem_info);
|
||||
ctx->forces.forces_mo_value_maxsize = mem_info.size;
|
||||
|
||||
if (forces_mo_value == NULL) {
|
||||
return qmckl_failwith( context,
|
||||
@ -5760,20 +5830,21 @@ integer(qmckl_exit_code) function qmckl_compute_forces_mo_value_doc(context, &
|
||||
|
||||
info = QMCKL_SUCCESS
|
||||
forces_mo_value = 0.0d0
|
||||
!do a=1,nucl_num
|
||||
! do j=1,point_num
|
||||
! do k=1,ao_num
|
||||
! c1 = forces_ao_value(k,j,1,a)
|
||||
! c2 = forces_ao_value(k,j,2,a)
|
||||
! c3 = forces_ao_value(k,j,3,a)
|
||||
! do i=1,mo_num
|
||||
! forces_mo_value(i,j,1,a) = forces_mo_value(i,j,1,a) + coefficient_t(i,k) * c1
|
||||
! forces_mo_value(i,j,2,a) = forces_mo_value(i,j,2,a) + coefficient_t(i,k) * c2
|
||||
! forces_mo_value(i,j,3,a) = forces_mo_value(i,j,3,a) + coefficient_t(i,k) * c3
|
||||
! end do
|
||||
! end do
|
||||
! end do
|
||||
!end do
|
||||
do a=1,nucl_num
|
||||
do j=1,point_num
|
||||
do k=1,ao_num
|
||||
c1 = forces_ao_value(k,j,1,a)
|
||||
c2 = forces_ao_value(k,j,2,a)
|
||||
c3 = forces_ao_value(k,j,3,a)
|
||||
if (c1 == 0.d0 .and. c2 == 0.d0 .and. c3 == 0.d0) cycle
|
||||
do i=1,mo_num
|
||||
forces_mo_value(i,j,1,a) = forces_mo_value(i,j,1,a) + coefficient_t(i,k) * c1
|
||||
forces_mo_value(i,j,2,a) = forces_mo_value(i,j,2,a) + coefficient_t(i,k) * c2
|
||||
forces_mo_value(i,j,3,a) = forces_mo_value(i,j,3,a) + coefficient_t(i,k) * c3
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
|
||||
!do a=1,nucl_num
|
||||
@ -5785,10 +5856,10 @@ integer(qmckl_exit_code) function qmckl_compute_forces_mo_value_doc(context, &
|
||||
! end do
|
||||
!end do
|
||||
|
||||
info = qmckl_dgemm(context,'N', 'N', mo_num, 3*nucl_num*point_num, ao_num, &
|
||||
1.0d0, coefficient_t(1,1), mo_num, &
|
||||
forces_ao_value(1,1,1,1), ao_num, &
|
||||
0.0d0, forces_mo_value(1,1,1,1), mo_num)
|
||||
! info = qmckl_dgemm(context,'N', 'N', mo_num, 3*nucl_num*point_num, ao_num, &
|
||||
! 1.0d0, coefficient_t(1,1), mo_num, &
|
||||
! forces_ao_value(1,1,1,1), ao_num, &
|
||||
! 0.0d0, forces_mo_value(1,1,1,1), mo_num)
|
||||
|
||||
|
||||
end function qmckl_compute_forces_mo_value_doc
|
||||
@ -6099,10 +6170,13 @@ integer function qmckl_compute_forces_mo_g_doc(context, &
|
||||
real (c_double ) , intent(in) :: ao_hessian(ao_num,4,point_num,3,nucl_num)
|
||||
real (c_double ) , intent(out) :: forces_mo_g(mo_num,3,point_num,3,nucl_num)
|
||||
integer*8 :: i,j,k, m, n,a
|
||||
double precision :: c1
|
||||
|
||||
|
||||
info = QMCKL_SUCCESS
|
||||
|
||||
|
||||
|
||||
forces_mo_g = 0.d0
|
||||
|
||||
!do j=1,point_num
|
||||
@ -6120,6 +6194,22 @@ integer function qmckl_compute_forces_mo_g_doc(context, &
|
||||
! end do
|
||||
!end do
|
||||
|
||||
do a=1,nucl_num
|
||||
do n = 1, 3
|
||||
do j=1,point_num
|
||||
do m = 1, 3
|
||||
do k=1,ao_num
|
||||
c1 = ao_hessian(k, m, j, n, a)
|
||||
if (c1 == 0.d0) cycle
|
||||
do i=1,mo_num
|
||||
forces_mo_g(i, m, j, n, a) = forces_mo_g(i, m, j, n, a) - coefficient_t(i,k) * c1
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
|
||||
|
||||
! do k=1,ao_num
|
||||
@ -6135,16 +6225,18 @@ integer function qmckl_compute_forces_mo_g_doc(context, &
|
||||
! ao_hessian(:, m, :, :, :), ao_num, &
|
||||
! 1.d0, forces_mo_g(:, m, :, :, :), mo_num)
|
||||
!end do
|
||||
do a=1,nucl_num
|
||||
do n = 1,3
|
||||
do m = 1, 3
|
||||
info = qmckl_dgemm(context,'N', 'N', mo_num, point_num, ao_num, &
|
||||
-1.d0, coefficient_t, mo_num, &
|
||||
ao_hessian(:, m, :, n, a), ao_num, &
|
||||
1.d0, forces_mo_g(:, m, :, n, a), mo_num)
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
|
||||
!do a=1,nucl_num
|
||||
! do n = 1,3
|
||||
! do m = 1, 3
|
||||
! info = qmckl_dgemm(context,'N', 'N', mo_num, point_num, ao_num, &
|
||||
! -1.d0, coefficient_t, mo_num, &
|
||||
! ao_hessian(:, m, :, n, a), ao_num, &
|
||||
! 1.d0, forces_mo_g(:, m, :, n, a), mo_num)
|
||||
! end do
|
||||
! end do
|
||||
!end do
|
||||
|
||||
|
||||
end function qmckl_compute_forces_mo_g_doc
|
||||
|
||||
File diff suppressed because it is too large
Load Diff
File diff suppressed because it is too large
Load Diff
@ -14,6 +14,7 @@ qmckl_determinant.org
|
||||
qmckl_sherman_morrison_woodbury.org
|
||||
qmckl_jastrow_champ.org
|
||||
qmckl_jastrow_champ_single.org
|
||||
qmckl_jastrow_champ_quad.org
|
||||
qmckl_local_energy.org
|
||||
qmckl_forces.org
|
||||
qmckl_utils.org
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user