From d36433cc1b9afa730c4e240449c60c974b49b78e Mon Sep 17 00:00:00 2001 From: Konstantin Karandashev Date: Sun, 28 Feb 2021 16:25:53 +0100 Subject: [PATCH 1/8] Introduced OpenMP parallelization and periodic boundary conditions support for FCHL19 routines. --- qml/kernels/fgradient_kernels.f90 | 6 +- qml/representations/facsf.f90 | 185 ++++++++++++++----------- qml/representations/representations.py | 19 ++- 3 files changed, 126 insertions(+), 84 deletions(-) diff --git a/qml/kernels/fgradient_kernels.f90 b/qml/kernels/fgradient_kernels.f90 index a65f8a57b..5bdb45f29 100644 --- a/qml/kernels/fgradient_kernels.f90 +++ b/qml/kernels/fgradient_kernels.f90 @@ -591,6 +591,7 @@ subroutine fatomic_local_gradient_kernel(x1, x2, dx2, q1, q2, n1, n2, nm1, nm2, sorted_derivs = 0.0d0 ! Presort the representation derivatives +!$OMP PARALLEL DO PRIVATE(idx2) schedule(dynamic) do b = 1, nm2 do i2 = 1, n2(b) idx2 = 0 @@ -606,10 +607,11 @@ subroutine fatomic_local_gradient_kernel(x1, x2, dx2, q1, q2, n1, n2, nm1, nm2, enddo enddo enddo +!$OMP END PARALLEL DO kernel = 0.0d0 - !$OMP PARALLEL DO PRIVATE(idx2_end,idx2_start,d,expd,idx1_start,idx1) schedule(dynamic) +!$OMP PARALLEL DO PRIVATE(idx2_end,idx2_start,d,expd,idx1_start,idx1) schedule(dynamic) do a = 1, nm1 idx1_start = sum(n1(:a)) - n1(a) + 1 @@ -642,7 +644,7 @@ subroutine fatomic_local_gradient_kernel(x1, x2, dx2, q1, q2, n1, n2, nm1, nm2, enddo enddo enddo - !$OMP END PARALLEL do +!$OMP END PARALLEL do deallocate(sorted_derivs) deallocate(d) diff --git a/qml/representations/facsf.f90 b/qml/representations/facsf.f90 index f4a8b5d81..25e587a98 100644 --- a/qml/representations/facsf.f90 +++ b/qml/representations/facsf.f90 @@ -634,7 +634,7 @@ end subroutine fgenerate_acsf_and_gradients subroutine fgenerate_fchl_acsf(coordinates, nuclear_charges, elements, & - & Rs2, Rs3, Ts, eta2, eta3, zeta, rcut, acut, natoms, rep_size, & + & Rs2, Rs3, Ts, eta2, eta3, zeta, rcut, acut, natoms, natoms_tot, rep_size, & & two_body_decay, three_body_decay, three_body_weight, rep) use acsf_utils, only: decay, calc_angle, calc_cos_angle @@ -652,20 +652,21 @@ subroutine fgenerate_fchl_acsf(coordinates, nuclear_charges, elements, & double precision, intent(in) :: zeta double precision, intent(in) :: rcut double precision, intent(in) :: acut - integer, intent(in) :: natoms + integer, intent(in) :: natoms, natoms_tot ! natoms_tot includes "virtual" atoms created to satisfy periodic boundary conditions. integer, intent(in) :: rep_size double precision, intent(in) :: two_body_decay double precision, intent(in) :: three_body_decay double precision, intent(in) :: three_body_weight double precision, intent(out), dimension(natoms, rep_size) :: rep +! Introduced to make OpenMP parallelization easier. + double precision, dimension(:, :, :), allocatable:: add_rep - integer :: i, j, k, l, n, m, o, p, q, s, z, nelements, nbasis2, nbasis3, nabasis + integer:: i, j, k, l, n, m, o, p, q, s, z, nelements, nbasis2, nbasis3, nabasis integer, allocatable, dimension(:) :: element_types double precision :: rij, rik, angle, cos_1, cos_2, cos_3, invcut - ! double precision :: angle_1, angle_2, angle_3 double precision, allocatable, dimension(:) :: radial, angular, a, b, c - double precision, allocatable, dimension(:, :) :: distance_matrix, rdecay, rep3 + double precision, allocatable, dimension(:, :) :: distance_matrix, rdecay double precision :: mu, sigma, ksi3 @@ -683,36 +684,36 @@ subroutine fgenerate_fchl_acsf(coordinates, nuclear_charges, elements, & ! number of element types nelements = size(elements) ! Allocate temporary - allocate(element_types(natoms)) + allocate(element_types(natoms_tot)) ! Store element index of every atom - ! !$OMP PARALLEL DO - do i = 1, natoms +!$OMP PARALLEL DO SCHEDULE(dynamic) + do i = 1, natoms_tot do j = 1, nelements - if (nuclear_charges(i) .eq. elements(j)) then + if (nuclear_charges(modulo(i-1, natoms)+1) .eq. elements(j)) then element_types(i) = j - continue + exit endif enddo enddo - ! !$OMP END PARALLEL DO +!$OMP END PARALLEL DO ! Get distance matrix ! Allocate temporary - allocate(distance_matrix(natoms, natoms)) + allocate(distance_matrix(natoms_tot, natoms_tot)) distance_matrix = 0.0d0 - ! !$OMP PARALLEL DO PRIVATE(rij) - do i = 1, natoms - do j = i+1, natoms +!$OMP PARALLEL DO PRIVATE(rij) SCHEDULE(dynamic) + do i = 1, natoms_tot + do j = i+1, natoms_tot rij = norm2(coordinates(j,:) - coordinates(i,:)) distance_matrix(i, j) = rij distance_matrix(j, i) = rij enddo enddo - ! !$OMP END PARALLEL DO +!$OMP END PARALLEL DO ! number of basis functions in the two body term nbasis2 = size(Rs2) @@ -721,17 +722,18 @@ subroutine fgenerate_fchl_acsf(coordinates, nuclear_charges, elements, & invcut = 1.0d0 / rcut ! pre-calculate the radial decay in the two body terms - rdecay = decay(distance_matrix, invcut, natoms) + rdecay = decay(distance_matrix, invcut, natoms_tot) ! Allocate temporary - allocate(radial(nbasis2)) + allocate(radial(nbasis2), add_rep(nbasis2, natoms, natoms)) radial = 0.0d0 - ! !$OMP PARALLEL DO PRIVATE(n,m,rij,radial) REDUCTION(+:rep) + add_rep=0.0d0 +!$OMP PARALLEL DO PRIVATE(n,m,rij,radial,mu,sigma) SCHEDULE(dynamic) do i = 1, natoms ! index of the element of atom i m = element_types(i) - do j = i + 1, natoms + do j = i + 1, natoms_tot ! index of the element of atom j n = element_types(j) ! distance between atoms i and j @@ -749,13 +751,23 @@ subroutine fgenerate_fchl_acsf(coordinates, nuclear_charges, elements, & enddo rep(i, (n-1)*nbasis2 + 1:n*nbasis2) = rep(i, (n-1)*nbasis2 + 1:n*nbasis2) + radial - rep(j, (m-1)*nbasis2 + 1:m*nbasis2) = rep(j, (m-1)*nbasis2 + 1:m*nbasis2) + radial +! rep(j, (m-1)*nbasis2 + 1:m*nbasis2) = rep(j, (m-1)*nbasis2 + 1:m*nbasis2) + radial + if (j<=natoms) add_rep(:, i, j)=radial endif enddo enddo - ! !$OMP END PARALLEL DO - +!$OMP END PARALLEL DO deallocate(radial) +!$OMP PARALLEL DO PRIVATE(m) SCHEDULE(dynamic) + do j = 2, natoms + do i = 1, j-1 + m = element_types(i) + rep(j, (m-1)*nbasis2 + 1:m*nbasis2) = rep(j, (m-1)*nbasis2 + 1:m*nbasis2) + add_rep(:, i, j) + enddo + enddo +!$OMP END PARALLEL DO + + deallocate(add_rep) ! number of radial basis functions in the three body term nbasis3 = size(Rs3) @@ -765,32 +777,26 @@ subroutine fgenerate_fchl_acsf(coordinates, nuclear_charges, elements, & ! Inverse of the three body cutoff invcut = 1.0d0 / acut ! pre-calculate the radial decay in the three body terms - rdecay = decay(distance_matrix, invcut, natoms) + rdecay = decay(distance_matrix, invcut, natoms_tot) ! Allocate temporary - allocate(rep3(natoms,rep_size)) allocate(a(3)) allocate(b(3)) allocate(c(3)) allocate(radial(nbasis3)) allocate(angular(nabasis)) - rep3 = 0.0d0 - - ! This could probably be done more efficiently if it's a bottleneck - ! Also the order is a bit wobbly compared to the tensorflow implementation - ! !$OMP PARALLEL DO PRIVATE(rij, n, rik, m, a, b, c, angle, radial, angular, & - ! !$OMP cos_1, cos_2, cos_3, mu, sigma, o, ksi3, & - ! !$OMP p, q, s, z) REDUCTION(+:rep3) COLLAPSE(2) SCHEDULE(dynamic) +!$OMP PARALLEL DO PRIVATE(rij, n, rik, m, a, b, c, angle, cos_1, cos_2, cos_3,& +!$OMP radial, ksi3, angular, o, p, q, s, z) SCHEDULE(dynamic) do i = 1, natoms - do j = 1, natoms - 1 + do j = 1, natoms_tot - 1 if (i .eq. j) cycle ! distance between atoms i and j rij = distance_matrix(i,j) if (rij > acut) cycle ! index of the element of atom j n = element_types(j) - do k = j + 1, natoms + do k = j + 1, natoms_tot if (i .eq. k) cycle if (j .eq. k) cycle ! distance between atoms i and k @@ -838,19 +844,16 @@ subroutine fgenerate_fchl_acsf(coordinates, nuclear_charges, elements, & ! calculate the indices that the three body terms should be added to z = s + (l-1) * nabasis ! Add the contributions from atoms i,j and k - rep3(i, z:z + nabasis - 1) = rep3(i, z:z + nabasis - 1) + angular * radial(l) * ksi3 + rep(i, z:z + nabasis - 1) = rep(i, z:z + nabasis - 1) + angular * radial(l) * ksi3 enddo enddo enddo enddo - ! !$OMP END PARALLEL DO - - rep = rep + rep3 +!$OMP END PARALLEL DO deallocate(element_types) deallocate(rdecay) deallocate(distance_matrix) - deallocate(rep3) deallocate(a) deallocate(b) deallocate(c) @@ -860,7 +863,7 @@ subroutine fgenerate_fchl_acsf(coordinates, nuclear_charges, elements, & end subroutine fgenerate_fchl_acsf subroutine fgenerate_fchl_acsf_and_gradients(coordinates, nuclear_charges, elements, & - & Rs2, Rs3, Ts, eta2, eta3, zeta, rcut, acut, natoms, rep_size, & + & Rs2, Rs3, Ts, eta2, eta3, zeta, rcut, acut, natoms, natoms_tot, rep_size, & & two_body_decay, three_body_decay, three_body_weight, rep, grad) use acsf_utils, only: decay, calc_angle, calc_cos_angle @@ -888,11 +891,15 @@ subroutine fgenerate_fchl_acsf_and_gradients(coordinates, nuclear_charges, eleme double precision, allocatable, dimension(:) :: exp_ln double precision, allocatable, dimension(:) :: log_Rs2 - integer, intent(in) :: natoms + integer, intent(in) :: natoms, natoms_tot ! natoms_tot includes "virtual" atoms created to satisfy periodic boundary conditions. integer, intent(in) :: rep_size double precision, intent(out), dimension(natoms, rep_size) :: rep double precision, intent(out), dimension(natoms, rep_size, natoms, 3) :: grad +! Introduced to make OpenMP parallelization easier. + double precision, dimension(:, :, :), allocatable:: add_rep + double precision, dimension(:, :, :, :), allocatable :: add_grad + integer :: i, j, k, l, m, n, p, q, s, t, z, nelements, nbasis2, nbasis3, nabasis, twobody_size integer, allocatable, dimension(:) :: element_types double precision :: rij, rik, angle, dot, rij2, rik2, invrij, invrik, invrij2, invrik2, invcut @@ -928,37 +935,37 @@ subroutine fgenerate_fchl_acsf_and_gradients(coordinates, nuclear_charges, eleme ! Number of unique elements nelements = size(elements) ! Allocate temporary - allocate(element_types(natoms)) + allocate(element_types(natoms_tot)) - ! Store element index of every atom - ! !$OMP PARALLEL DO - do i = 1, natoms +! Store element index of every atom +!$OMP PARALLEL DO SCHEDULE(dynamic) + do i = 1, natoms_tot do j = 1, nelements - if (nuclear_charges(i) .eq. elements(j)) then + if (nuclear_charges(modulo(i-1, natoms)+1) .eq. elements(j)) then element_types(i) = j - continue + exit endif enddo enddo - ! !$OMP END PARALLEL DO +!$OMP END PARALLEL DO ! Get distance matrix ! Allocate temporary - allocate(distance_matrix(natoms, natoms)) - allocate(sq_distance_matrix(natoms, natoms)) - allocate(inv_distance_matrix(natoms, natoms)) - allocate(inv_sq_distance_matrix(natoms, natoms)) + allocate(distance_matrix(natoms_tot, natoms_tot)) + allocate(sq_distance_matrix(natoms_tot, natoms_tot)) + allocate(inv_distance_matrix(natoms_tot, natoms_tot)) + allocate(inv_sq_distance_matrix(natoms_tot, natoms_tot)) distance_matrix = 0.0d0 sq_distance_matrix = 0.0d0 inv_distance_matrix = 0.0d0 inv_sq_distance_matrix = 0.0d0 - ! !$OMP PARALLEL DO PRIVATE(rij,rij2,invrij,invrij2) SCHEDULE(dynamic) - do i = 1, natoms - do j = i+1, natoms +!$OMP PARALLEL DO PRIVATE(rij,rij2,invrij,invrij2) SCHEDULE(dynamic) + do i = 1, natoms_tot + do j = i+1, natoms_tot rij = norm2(coordinates(j,:) - coordinates(i,:)) distance_matrix(i, j) = rij distance_matrix(j, i) = rij @@ -973,7 +980,7 @@ subroutine fgenerate_fchl_acsf_and_gradients(coordinates, nuclear_charges, eleme inv_sq_distance_matrix(j, i) = invrij2 enddo enddo - ! !$OMP END PARALLEL DO +!$OMP END PARALLEL DO ! Number of two body basis functions @@ -982,7 +989,7 @@ subroutine fgenerate_fchl_acsf_and_gradients(coordinates, nuclear_charges, eleme ! Inverse of the two body cutoff distance invcut = 1.0d0 / rcut ! Pre-calculate the two body decay - rdecay = decay(distance_matrix, invcut, natoms) + rdecay = decay(distance_matrix, invcut, natoms_tot) ! Allocate temporary allocate(radial_base(nbasis2)) @@ -995,14 +1002,17 @@ subroutine fgenerate_fchl_acsf_and_gradients(coordinates, nuclear_charges, eleme grad =0.0d0 rep = 0.0d0 + allocate(add_rep(nbasis2, natoms, natoms), add_grad(nbasis2, 3, natoms, natoms)) + add_rep=0.0d0 + add_grad=0.0d0 log_Rs2(:) = log(Rs2(:)) - ! !$OMP PARALLEL DO PRIVATE(m,n,rij,invrij,radial_base,radial,radial_part,part) REDUCTION(+:rep,grad) & - ! !$OMP SCHEDULE(dynamic) +!$OMP PARALLEL DO PRIVATE(m, n, rij, invrij, mu, sigma, exp_s2, exp_ln,& +!$OMP scaling, radial_base, radial, dx, part, dscal, ddecay) SCHEDULE(dynamic) do i = 1, natoms ! The element index of atom i m = element_types(i) - do j = i + 1, natoms + do j = i + 1, natoms_tot ! The element index of atom j n = element_types(j) ! Distance between atoms i and j @@ -1023,8 +1033,8 @@ subroutine fgenerate_fchl_acsf_and_gradients(coordinates, nuclear_charges, eleme radial(:) = radial_base(:) * scaling * rdecay(i,j) rep(i, (n-1)*nbasis2 + 1:n*nbasis2) = rep(i, (n-1)*nbasis2 + 1:n*nbasis2) + radial - rep(j, (m-1)*nbasis2 + 1:m*nbasis2) = rep(j, (m-1)*nbasis2 + 1:m*nbasis2) + radial - +! rep(j, (m-1)*nbasis2 + 1:m*nbasis2) = rep(j, (m-1)*nbasis2 + 1:m*nbasis2) + radial + if (j<=natoms) add_rep(:, i, j)=radial do k = 1, 3 @@ -1043,15 +1053,31 @@ subroutine fgenerate_fchl_acsf_and_gradients(coordinates, nuclear_charges, eleme ! The gradients wrt coordinates grad(i, (n-1)*nbasis2 + 1:n*nbasis2, i, k) = grad(i, (n-1)*nbasis2 + 1:n*nbasis2, i, k) + part - grad(i, (n-1)*nbasis2 + 1:n*nbasis2, j, k) = grad(i, (n-1)*nbasis2 + 1:n*nbasis2, j, k) - part - grad(j, (m-1)*nbasis2 + 1:m*nbasis2, j, k) = grad(j, (m-1)*nbasis2 + 1:m*nbasis2, j, k) - part - grad(j, (m-1)*nbasis2 + 1:m*nbasis2, i, k) = grad(j, (m-1)*nbasis2 + 1:m*nbasis2, i, k) + part - + if (j<=natoms) then + grad(i, (n-1)*nbasis2 + 1:n*nbasis2, j, k) = grad(i, (n-1)*nbasis2 + 1:n*nbasis2, j, k) - part + grad(j, (m-1)*nbasis2 + 1:m*nbasis2, i, k) = grad(j, (m-1)*nbasis2 + 1:m*nbasis2, i, k) + part +! grad(j, (m-1)*nbasis2 + 1:m*nbasis2, j, k) = grad(j, (m-1)*nbasis2 + 1:m*nbasis2, j, k) - part + add_grad(:, k, i, j)=part + endif enddo endif enddo enddo - ! !$OMP END PARALLEL DO +!$OMP END PARALLEL DO + +!$OMP PARALLEL DO PRIVATE(m) SCHEDULE(dynamic) + do j = 2, natoms + do i = 1, j-1 + m = element_types(i) + rep(j, (m-1)*nbasis2 + 1:m*nbasis2) = rep(j, (m-1)*nbasis2 + 1:m*nbasis2) + add_rep(:, i, j) + do k=1, 3 + grad(j, (m-1)*nbasis2 + 1:m*nbasis2, j, k) = grad(j, (m-1)*nbasis2 + 1:m*nbasis2, j, k) - add_grad(:, k, i, j) + enddo + enddo + enddo +!$OMP END PARALLEL DO + + deallocate(add_rep, add_grad) deallocate(radial_base) deallocate(radial) @@ -1069,7 +1095,7 @@ subroutine fgenerate_fchl_acsf_and_gradients(coordinates, nuclear_charges, eleme ! Inverse of the three body cutoff distance invcut = 1.0d0 / acut ! Pre-calculate the three body decay - rdecay = decay(distance_matrix, invcut, natoms) + rdecay = decay(distance_matrix, invcut, natoms_tot) ! Allocate temporary allocate(atom_rep(rep_size - twobody_size)) @@ -1114,17 +1140,17 @@ subroutine fgenerate_fchl_acsf_and_gradients(coordinates, nuclear_charges, eleme allocate(d_atm_extra_j(3)) allocate(d_atm_extra_k(3)) - ! ! This could probably be done more efficiently if it's a bottleneck - ! ! The order is a bit wobbly compared to the tensorflow implementation - ! !$OMP PARALLEL DO PRIVATE(atom_rep,atom_grad,a,b,c,radial,angular_base, & - ! !$OMP angular,d_angular,rij,n,rij2,invrij,invrij2,d_angular_d_i, & - ! !$OMP d_angular_d_j,d_angular_d_k,rik,m,rik2,invrik,invrik2,angle, & - ! !$OMP p,q,dot,d_radial,d_radial_d_i,d_radial_d_j,d_radial_d_k,s,z, & - ! !$OMP d_ijdecay,d_ikdecay) SCHEDULE(dynamic) +!$OMP PARALLEL DO PRIVATE(atom_rep, atom_grad, rij, n, rij2, invrij, invrij2,& +!$OMP rik, m, rik2, invrik, invrjk, invrik2, a, b, c, angle, cos_i, cos_k,& +!$OMP cos_j, radial_base, radial, p, q, dot, angular, d_angular, d_angular_d_j,& +!$OMP d_angular_d_k, d_angular_d_i, d_radial, d_radial_d_j, d_radial_d_k,& +!$OMP d_radial_d_i, d_ijdecay, d_ikdecay, invr_atm, atm, atm_i, atm_j, atm_k, vi,& +!$OMP vj, vk, d_atm_ii, d_atm_ij, d_atm_ik, d_atm_ji, d_atm_jj, d_atm_jk, d_atm_ki,& +!$OMP d_atm_kj, d_atm_kk, d_atm_extra_i, d_atm_extra_j, d_atm_extra_k, s, z) SCHEDULE(dynamic) do i = 1, natoms atom_rep = 0.0d0 atom_grad = 0.0d0 - do j = 1, natoms - 1 + do j = 1, natoms_tot - 1 if (i .eq. j) cycle ! distance between atom i and j rij = distance_matrix(i,j) @@ -1137,7 +1163,7 @@ subroutine fgenerate_fchl_acsf_and_gradients(coordinates, nuclear_charges, eleme invrij = inv_distance_matrix(i,j) ! inverse squared distance between atom i and j invrij2 = inv_sq_distance_matrix(i,j) - do k = j + 1, natoms + do k = j + 1, natoms_tot if (i .eq. k) cycle ! distance between atom i and k rik = distance_matrix(i,k) @@ -1250,8 +1276,8 @@ subroutine fgenerate_fchl_acsf_and_gradients(coordinates, nuclear_charges, eleme & angular * radial(l) * (atm_i * d_atm_ii(t) + atm_j * d_atm_ij(t) & & + atm_k * d_atm_ik(t) + d_atm_extra_i(t)) * three_body_weight * rdecay(i,j) * rdecay(i,k) + & & angular * radial(l) * (d_ijdecay(t) * rdecay(i,k) + rdecay(i,j) * d_ikdecay(t)) * atm - ! Add up all gradient contributions wrt atom j + if (j<=natoms) & atom_grad(z:z + nabasis - 1, j, t) = atom_grad(z:z + nabasis - 1, j, t) + & & d_angular * d_angular_d_j(t) * radial(l) * atm * rdecay(i,j) * rdecay(i,k) + & & angular * d_radial(l) * d_radial_d_j(t) * atm * rdecay(i,j) * rdecay(i,k) + & @@ -1260,6 +1286,7 @@ subroutine fgenerate_fchl_acsf_and_gradients(coordinates, nuclear_charges, eleme & angular * radial(l) * d_ijdecay(t) * rdecay(i,k) * atm ! Add up all gradient contributions wrt atom k + if (k<=natoms) & atom_grad(z:z + nabasis - 1, k, t) = atom_grad(z:z + nabasis - 1, k, t) + & & d_angular * d_angular_d_k(t) * radial(l) * atm * rdecay(i,j) * rdecay(i,k) + & & angular * d_radial(l) * d_radial_d_k(t) * atm * rdecay(i,j) * rdecay(i,k) + & @@ -1274,7 +1301,7 @@ subroutine fgenerate_fchl_acsf_and_gradients(coordinates, nuclear_charges, eleme rep(i, twobody_size + 1:) = rep(i, twobody_size + 1:) + atom_rep grad(i, twobody_size + 1:,:,:) = grad(i, twobody_size + 1:,:,:) + atom_grad enddo - ! !$OMP END PARALLEL DO +!$OMP END PARALLEL DO deallocate(rdecay) deallocate(element_types) diff --git a/qml/representations/representations.py b/qml/representations/representations.py index 9cf89c360..a8f2b5386 100644 --- a/qml/representations/representations.py +++ b/qml/representations/representations.py @@ -630,7 +630,7 @@ def generate_acsf(nuclear_charges, coordinates, elements = [1,6,7,8,16], nRs2 = def generate_fchl_acsf(nuclear_charges, coordinates, elements = [1,6,7,8,16], nRs2=24, nRs3=20, nFourier=1, eta2=0.32, eta3=2.7, zeta=np.pi, rcut=8.0, acut=8.0, two_body_decay=1.8, three_body_decay=0.57, three_body_weight=13.4, - pad=False, gradients=False): + pad=False, gradients=False, cell=None): """ FCHL-ACSF @@ -667,6 +667,8 @@ def generate_fchl_acsf(nuclear_charges, coordinates, elements = [1,6,7,8,16], :type gradients: boolean :return: Atom-centered symmetry functions representation :rtype: numpy array + :param cell: Unit cell vectors. The presence of this keyword argument will generate a periodic representation. + :type cell: numpy array """ Rs2 = np.linspace(0, rcut, 1+nRs2)[1:] @@ -681,11 +683,22 @@ def generate_fchl_acsf(nuclear_charges, coordinates, elements = [1,6,7,8,16], # Normalization constant for three-body three_body_weight = np.sqrt(eta3/np.pi) * three_body_weight + natoms_tot=natoms + if cell is not None: + nExtend = (np.floor(rcut/np.linalg.norm(cell,2,axis = 0)) + 1).astype(int) + true_coords=coordinates + for i in range(-nExtend[0],nExtend[0] + 1): + for j in range(-nExtend[1],nExtend[1] + 1): + for k in range(-nExtend[2],nExtend[2] + 1): + if not (i == 0 and j == 0 and k == 0): + true_coords = np.append(true_coords, coordinates + i*cell[0,:] + j*cell[1,:] + k*cell[2,:], axis=0) + natoms_tot+=natoms + coordinates=true_coords if gradients is False: rep = fgenerate_fchl_acsf(coordinates, nuclear_charges, elements, Rs2, Rs3, - Ts, eta2, eta3, zeta, rcut, acut, natoms, descr_size, + Ts, eta2, eta3, zeta, rcut, acut, natoms, natoms_tot, descr_size, two_body_decay, three_body_decay, three_body_weight) if pad is not False: @@ -705,7 +718,7 @@ def generate_fchl_acsf(nuclear_charges, coordinates, elements = [1,6,7,8,16], exit() (rep, grad) = fgenerate_fchl_acsf_and_gradients(coordinates, nuclear_charges, elements, Rs2, Rs3, - Ts, eta2, eta3, zeta, rcut, acut, natoms, descr_size, + Ts, eta2, eta3, zeta, rcut, acut, natoms, natoms_tot, descr_size, two_body_decay, three_body_decay, three_body_weight) From e40699bca152aaad6776480170aa620783ddc1fb Mon Sep 17 00:00:00 2001 From: Max Schwilk Date: Thu, 4 Mar 2021 17:01:55 +0100 Subject: [PATCH 2/8] first version of test job for PBC with FCHL19 forces learning --- test/data/LiH_crystal_all.xyz | 152 +++++++++++++++++++++++++++++ test/data/LiH_crystal_test.xyz | 32 ++++++ test/data/LiH_crystal_training.xyz | 40 ++++++++ test/test_fchl_acsf_forces.py | 96 +++++++++++++++++- 4 files changed, 319 insertions(+), 1 deletion(-) create mode 100644 test/data/LiH_crystal_all.xyz create mode 100644 test/data/LiH_crystal_test.xyz create mode 100644 test/data/LiH_crystal_training.xyz diff --git a/test/data/LiH_crystal_all.xyz b/test/data/LiH_crystal_all.xyz new file mode 100644 index 000000000..aa888de73 --- /dev/null +++ b/test/data/LiH_crystal_all.xyz @@ -0,0 +1,152 @@ +2 +Properties=species:S:1:pos:R:3:forces:R:3 Lattice="2.0 0.0 0.0 0.0 2.0 0.0 0.0 0.0 2.0" energy=-210.34143869968557 +H 1.7469600972641757 1.3105406236661603 1.8078885055514737 -0.05515010202844178 -0.016138638841873887 -0.07494808520830193 +Li 1.3428892907885033 0.34725373984758257 1.1993919048712993 0.05608062426846 0.019533274244299792 0.07536758225127621 +2 +Properties=species:S:1:pos:R:3:forces:R:3 Lattice="2.0 0.0 0.0 0.0 2.0 0.0 0.0 0.0 2.0" energy=-208.1561426455114 +H 0.5421120831942108 1.123751314386705 1.850047836343 0.28176960731737 -0.04513199169175898 0.25950818032265166 +Li 1.228635938071446 1.032454202454703 0.3783686165462776 -0.28071508108709975 0.045080998695914434 -0.26245776947660476 +2 +Properties=species:S:1:pos:R:3:forces:R:3 Lattice="2.0 0.0 0.0 0.0 2.0 0.0 0.0 0.0 2.0" energy=-210.09645903521653 +H 0.027158196041101146 1.6731862097711294 1.6134366515700422 0.10292634370603265 0.1275333232368676 0.04737263331395791 +Li 0.7941335004585335 0.3699666635090202 1.8523215817301748 -0.10497046544138594 -0.12381993344976971 -0.04826856097448054 +2 +Properties=species:S:1:pos:R:3:forces:R:3 Lattice="2.0 0.0 0.0 0.0 2.0 0.0 0.0 0.0 2.0" energy=-205.90144564737625 +H 1.8388327791612642 0.9993481253471506 0.4385729436404704 -0.22101994991201718 -0.16858417104422457 0.34610561623965597 +Li 1.4674006284615426 0.6984981715831435 1.1253342276681788 0.22287617184884612 0.17892796089444607 -0.34832478273018663 +2 +Properties=species:S:1:pos:R:3:forces:R:3 Lattice="2.0 0.0 0.0 0.0 2.0 0.0 0.0 0.0 2.0" energy=-210.78178622831706 +H 1.0168310046613556 0.8612585238423034 0.0009601399098819741 -0.04918730902954849 -0.0494842923887146 -0.002704796073556892 +Li 0.44419313332873056 0.2832776533608796 0.9973760638707188 0.05166094970816282 0.05193243921592955 -0.0008919259574161065 +2 +Properties=species:S:1:pos:R:3:forces:R:3 Lattice="2.0 0.0 0.0 0.0 2.0 0.0 0.0 0.0 2.0" energy=-204.93288634154862 +H 1.6369978211256002 0.3135443027001339 0.17708217034721807 0.286981721245257 0.5591418497105554 -0.3393644337261212 +Li 1.9117219165709483 0.8863594602481111 1.8673689857714335 -0.2855514427745901 -0.5622380627972693 0.33926335540009644 +2 +Properties=species:S:1:pos:R:3:forces:R:3 Lattice="2.0 0.0 0.0 0.0 2.0 0.0 0.0 0.0 2.0" energy=-210.73524383643314 +H 1.7700764233531123 0.7893408914492881 1.088955393401397 -0.048281281265649345 0.05406751470004684 0.0033540969189816416 +Li 1.2528169082220253 1.403038973917562 0.09006815480692909 0.05071086876279346 -0.05608528360044185 0.00022152166059696832 +2 +Properties=species:S:1:pos:R:3:forces:R:3 Lattice="2.0 0.0 0.0 0.0 2.0 0.0 0.0 0.0 2.0" energy=-210.7055029075143 +H 1.6106685007271484 1.0136065810407606 1.3399205019703946 0.07442600324468429 -0.0551332730156715 0.06893938171309927 +Li 0.22790861256500783 0.45451125078373633 0.13344906283841884 -0.07884652763782285 0.05491810155697552 -0.07288355009254194 +2 +Properties=species:S:1:pos:R:3:forces:R:3 Lattice="2.0 0.0 0.0 0.0 2.0 0.0 0.0 0.0 2.0" energy=-211.46640045705647 +H 1.8871080860691987 0.8752458196713906 0.6103866813399352 0.027083053922416 0.0061748501139223255 0.0009190170012606802 +Li 0.7841608913523739 0.008999028007995014 1.3435906437324703 -0.022815401313884376 0.0020283655074712637 -0.009497866107779085 +2 +Properties=species:S:1:pos:R:3:forces:R:3 Lattice="2.0 0.0 0.0 0.0 2.0 0.0 0.0 0.0 2.0" energy=-208.31080335905835 +H 0.751218546267328 1.8293859346784571 0.8411942880456913 0.10733713518674659 -0.1813798156351406 -0.18119888599741052 +Li 1.1181411223414048 1.2045494112405455 0.21726908998803207 -0.11639787090312986 0.18193530292844318 0.18178408033215054 +2 +Properties=species:S:1:pos:R:3:forces:R:3 Lattice="2.0 0.0 0.0 0.0 2.0 0.0 0.0 0.0 2.0" energy=-209.77427812121766 +H 1.5755046749473112 1.4508766716642294 0.9255583812921313 -0.10274734710384964 -0.032219684225297446 0.04469996473204668 +Li 0.9410022181884197 1.288271197130342 1.8422141803657355 0.10418706186491483 0.03395014705076274 -0.04775642087221893 +2 +Properties=species:S:1:pos:R:3:forces:R:3 Lattice="2.0 0.0 0.0 0.0 2.0 0.0 0.0 0.0 2.0" energy=-207.4907749598616 +H 1.9765799421188173 0.2949641306899158 0.6665239835320342 0.11701323176693335 0.2640164206981671 0.12532456152730848 +Li 0.14861512557309342 1.1299792755768918 0.8832177768619527 -0.1289002629261331 -0.267848137090989 -0.1230727462782803 +2 +Properties=species:S:1:pos:R:3:forces:R:3 Lattice="2.0 0.0 0.0 0.0 2.0 0.0 0.0 0.0 2.0" energy=-203.76022483437532 +H 1.7832481213684168 1.1384329182262618 0.028729360472953713 0.013126503911447061 -0.5845541870211601 -0.6338906563014737 +Li 1.7925113387749048 0.6991376801580882 1.537288146053771 -0.012355239328728135 0.5808077238808345 0.6409913147853306 +2 +Properties=species:S:1:pos:R:3:forces:R:3 Lattice="2.0 0.0 0.0 0.0 2.0 0.0 0.0 0.0 2.0" energy=-209.62779094440157 +H 0.7610588112475634 0.24461924434332794 1.634205810533692 -0.11273891948004727 -0.11632670596923639 0.13849374657889763 +Li 0.37212249222327 1.880627601969749 0.48462734990577383 0.10772952373292766 0.12424046313712078 -0.13527630474529206 +2 +Properties=species:S:1:pos:R:3:forces:R:3 Lattice="2.0 0.0 0.0 0.0 2.0 0.0 0.0 0.0 2.0" energy=-210.47076717382572 +H 0.2273823110242048 0.5559492502404069 1.637836758506292 -0.0851517715922081 -0.11531125765937161 -0.01426116804189391 +Li 1.3591633348032384 1.8662612430851568 1.5317806392640836 0.07879869554202434 0.11175371638784337 0.01872839925449693 +2 +Properties=species:S:1:pos:R:3:forces:R:3 Lattice="2.0 0.0 0.0 0.0 2.0 0.0 0.0 0.0 2.0" energy=-210.33446908697192 +H 0.5169590753697664 0.7416018162127209 0.28258963145894955 -0.042598862667708126 0.06979988236433743 0.038964187468775446 +Li 0.17483139755064836 1.4429162784323317 1.1675869685949671 0.04469666886338697 -0.07150518510612097 -0.038671505395925954 +2 +Properties=species:S:1:pos:R:3:forces:R:3 Lattice="2.0 0.0 0.0 0.0 2.0 0.0 0.0 0.0 2.0" energy=-205.11099128069407 +H 1.546344760224777 0.4602754563883884 1.833051407765566 -0.1596347733516269 0.4722982331388814 0.527863271499883 +Li 1.4006076540247496 0.912105388620791 0.35829211172807973 0.15562273655359407 -0.4679436737919508 -0.5314842153952676 +2 +Properties=species:S:1:pos:R:3:forces:R:3 Lattice="2.0 0.0 0.0 0.0 2.0 0.0 0.0 0.0 2.0" energy=-210.19627671692353 +H 1.2488266976231195 1.781675227044553 0.048209683160555405 0.07225356081485984 0.06302750536400226 0.009499287041346172 +Li 1.851533708686724 0.13039164075443477 1.0248588282589302 -0.07332310988078239 -0.07104735860181444 -0.0090847059079322 +2 +Properties=species:S:1:pos:R:3:forces:R:3 Lattice="2.0 0.0 0.0 0.0 2.0 0.0 0.0 0.0 2.0" energy=-208.22815198343892 +H 0.6660709328369625 0.6111848187354765 1.8993235598624758 -0.20092080259819828 0.24603601627548255 0.1135129361662014 +Li 0.15484187038923536 1.3163448968315028 0.13380429858981313 0.20255665461548672 -0.24846209057675384 -0.1255900417451442 +2 +Properties=species:S:1:pos:R:3:forces:R:3 Lattice="2.0 0.0 0.0 0.0 2.0 0.0 0.0 0.0 2.0" energy=-207.58115625484743 +H 0.33520134576418403 0.6498328960624737 1.97595666775638 0.2683683904298222 -0.4231543730218725 0.04005770247610885 +Li 0.7090217441607694 1.9479248003510565 0.010424693549549335 -0.2756287029445859 0.4229944309513502 -0.05541006305092011 +2 +Properties=species:S:1:pos:R:3:forces:R:3 Lattice="2.0 0.0 0.0 0.0 2.0 0.0 0.0 0.0 2.0" energy=-205.68766305590884 +H 0.6991987735247565 1.6037349788530562 1.2208702649325034 0.19473006228560652 -0.09445176695463661 0.3550448721773445 +Li 1.0041615376028266 1.4550517768499116 1.9892517660494589 -0.20079496138786623 0.10046342012006015 -0.35492625682796397 +2 +Properties=species:S:1:pos:R:3:forces:R:3 Lattice="2.0 0.0 0.0 0.0 2.0 0.0 0.0 0.0 2.0" energy=-146.8630693271275 +H 0.3103491102517437 1.241938223056992 1.072729408425892 -7.916938547982216 1.145367280941611 -0.0568068823103567 +Li 0.04834006988435835 1.2796657700568523 1.0708485433222104 7.907393602296153 -1.1585550666420672 0.057028535670367544 +2 +Properties=species:S:1:pos:R:3:forces:R:3 Lattice="2.0 0.0 0.0 0.0 2.0 0.0 0.0 0.0 2.0" energy=-210.84937993911183 +H 0.9826396567525186 0.3956944942235976 1.0293187007718965 0.007737459785600753 0.045633173694166806 0.03540358077018818 +Li 1.950682311583713 1.1321473282639127 1.4603857837198702 -0.010309551820136265 -0.048826449253047466 -0.03892603032439568 +2 +Properties=species:S:1:pos:R:3:forces:R:3 Lattice="2.0 0.0 0.0 0.0 2.0 0.0 0.0 0.0 2.0" energy=-208.41393938607132 +H 1.611810633906867 0.2519884930657792 0.5134967256789191 -0.2445149123942688 -0.21597810131521233 0.08123472099401788 +Li 0.8652482863697526 1.7880942947688458 0.702790881590545 0.24401220796303003 0.21973869068084012 -0.08287445706085114 +2 +Properties=species:S:1:pos:R:3:forces:R:3 Lattice="2.0 0.0 0.0 0.0 2.0 0.0 0.0 0.0 2.0" energy=-211.48912509225772 +H 0.48757626410196586 0.5611662537358879 1.9635928646173995 -0.026789830059731734 -0.0047521346621653415 0.0056915641190005695 +Li 1.5920828553293533 1.3695590558128652 1.1381109198893085 0.022390216998396273 -0.004036860331450565 0.0030402851052320212 +2 +Properties=species:S:1:pos:R:3:forces:R:3 Lattice="2.0 0.0 0.0 0.0 2.0 0.0 0.0 0.0 2.0" energy=-205.402618904674 +H 0.08705706944387637 1.9541357163145763 0.8148046208845077 0.15374240109363915 -0.006771054467366847 0.3549970889851254 +Li 0.30894099679302744 1.9446523651335261 1.624224289234903 -0.1622405655933489 0.0072524370011842025 -0.35418239662966344 +2 +Properties=species:S:1:pos:R:3:forces:R:3 Lattice="2.0 0.0 0.0 0.0 2.0 0.0 0.0 0.0 2.0" energy=-210.88050011212317 +H 0.6091167182824828 1.3611434989785527 0.6970211853178458 -0.04963684608722568 0.05023404972749873 -0.048415546988366864 +Li 1.7251811564142063 0.018189163332650304 1.9703843514414876 0.033865825316945375 -0.04403280355970984 0.042636819619115474 +2 +Properties=species:S:1:pos:R:3:forces:R:3 Lattice="2.0 0.0 0.0 0.0 2.0 0.0 0.0 0.0 2.0" energy=-210.77144264277774 +H 0.8932071011537388 1.2981984032516969 1.6939764959608663 0.037269032997885254 0.03772618762294472 0.08358864533807586 +Li 1.238075593194045 0.23854138005793213 0.38774468705418297 -0.04032171198505258 -0.03142958780248023 -0.08121544250454099 +2 +Properties=species:S:1:pos:R:3:forces:R:3 Lattice="2.0 0.0 0.0 0.0 2.0 0.0 0.0 0.0 2.0" energy=-209.52240451811292 +H 1.2081248409776497 1.8763329475126198 0.5245989961237525 0.13308330732558543 0.1357256725127428 0.1208069585595496 +Li 1.978741352030817 0.2887620867141434 1.0321698160632717 -0.13148742364419697 -0.12080999689923766 -0.12237640597062588 +2 +Properties=species:S:1:pos:R:3:forces:R:3 Lattice="2.0 0.0 0.0 0.0 2.0 0.0 0.0 0.0 2.0" energy=-207.07992135496852 +H 1.4313121083301743 0.167230811955142 1.0654810569217192 0.3845911913724531 0.12699191634917872 -0.20234159781489036 +Li 0.18026189432537443 0.36509977389329573 0.784417332212211 -0.3808170960555125 -0.1276519344086368 0.19244070839148275 +2 +Properties=species:S:1:pos:R:3:forces:R:3 Lattice="2.0 0.0 0.0 0.0 2.0 0.0 0.0 0.0 2.0" energy=-209.42471109764404 +H 0.2863550557467962 0.5403983677631059 1.2091050222686635 -0.17938746758451524 -0.13594357769769283 -0.1518444643231518 +Li 1.7013527171404086 0.037827601853055004 0.6271570433726206 0.17688963920673656 0.13595493733515807 0.15079472377258263 +2 +Properties=species:S:1:pos:R:3:forces:R:3 Lattice="2.0 0.0 0.0 0.0 2.0 0.0 0.0 0.0 2.0" energy=-210.8464400606922 +H 1.115631077279221 0.6107122764944715 0.11832655553829352 -0.020643087175397368 0.025543691133341373 -0.044400612691328445 +Li 0.9154628943603618 1.4920148173736698 1.24246582256913 0.02194696422856776 -0.03306740536831865 0.040629256614064024 +2 +Properties=species:S:1:pos:R:3:forces:R:3 Lattice="2.0 0.0 0.0 0.0 2.0 0.0 0.0 0.0 2.0" energy=-210.65409362767662 +H 1.958764583801193 1.618254016886918 1.67088147423364 0.04710863957682776 -0.03441263562899041 0.09457839337754242 +Li 0.8890679925379519 1.3501094012154145 0.356523142900558 -0.04007135838707443 0.03718283478467768 -0.09125227887537685 +2 +Properties=species:S:1:pos:R:3:forces:R:3 Lattice="2.0 0.0 0.0 0.0 2.0 0.0 0.0 0.0 2.0" energy=-198.27023740453092 +H 1.2841025399739825 1.0560738997199883 0.4076419004771905 -0.6622270842275676 -0.6621719165769338 0.29120726196817465 +Li 0.862900092659493 0.6349124368264709 0.5833565026928345 0.6623678945012492 0.6623120120756423 -0.2882727137301109 +2 +Properties=species:S:1:pos:R:3:forces:R:3 Lattice="2.0 0.0 0.0 0.0 2.0 0.0 0.0 0.0 2.0" energy=-210.62932122368318 +H 1.9990647323447277 0.953558635929386 0.7911340788025603 -0.05594055535694503 -0.057208550417013004 -0.05611647378515233 +Li 1.3290457311009285 0.2507507741141308 0.11711990229764058 0.059840751513080936 0.06045633930206612 0.05994014971060668 +2 +Properties=species:S:1:pos:R:3:forces:R:3 Lattice="2.0 0.0 0.0 0.0 2.0 0.0 0.0 0.0 2.0" energy=-211.11569239713003 +H 0.3293052666603564 0.6745606943807767 1.8083121053115703 0.020429525340863397 -0.02662672320307269 0.005988971418682931 +Li 0.6514315935216397 1.7481427802214673 0.8452546116538187 -0.02444469103096769 0.023429251001441775 0.0019257157545129466 +2 +Properties=species:S:1:pos:R:3:forces:R:3 Lattice="2.0 0.0 0.0 0.0 2.0 0.0 0.0 0.0 2.0" energy=-207.53626571234986 +H 0.8566887767663101 0.514330032046471 0.9820183385163381 -0.23398535394816644 0.194856087268374 0.07061133642254105 +Li 0.10276067711712411 1.0092298537393227 1.167843352890332 0.233035325671203 -0.19619665016142407 -0.07817745982712898 +2 +Properties=species:S:1:pos:R:3:forces:R:3 Lattice="2.0 0.0 0.0 0.0 2.0 0.0 0.0 0.0 2.0" energy=-206.9335273563688 +H 0.9365375039915607 1.070921326767058 1.066203800139247 0.2107657636447698 -0.10299811558112562 0.2770384753981818 +Li 1.3880639850834262 0.836584073664475 1.791586511019425 -0.21262545187531945 0.11296453822579205 -0.27676809894460996 diff --git a/test/data/LiH_crystal_test.xyz b/test/data/LiH_crystal_test.xyz new file mode 100644 index 000000000..52c9facbe --- /dev/null +++ b/test/data/LiH_crystal_test.xyz @@ -0,0 +1,32 @@ +2 +Properties=species:S:1:pos:R:3:forces:R:3 Lattice="2.0 0.0 0.0 0.0 2.0 0.0 0.0 0.0 2.0" energy=-209.42471109764404 +H 0.2863550557467962 0.5403983677631059 1.2091050222686635 -0.17938746758451524 -0.13594357769769283 -0.1518444643231518 +Li 1.7013527171404086 0.037827601853055004 0.6271570433726206 0.17688963920673656 0.13595493733515807 0.15079472377258263 +2 +Properties=species:S:1:pos:R:3:forces:R:3 Lattice="2.0 0.0 0.0 0.0 2.0 0.0 0.0 0.0 2.0" energy=-210.8464400606922 +H 1.115631077279221 0.6107122764944715 0.11832655553829352 -0.020643087175397368 0.025543691133341373 -0.044400612691328445 +Li 0.9154628943603618 1.4920148173736698 1.24246582256913 0.02194696422856776 -0.03306740536831865 0.040629256614064024 +2 +Properties=species:S:1:pos:R:3:forces:R:3 Lattice="2.0 0.0 0.0 0.0 2.0 0.0 0.0 0.0 2.0" energy=-210.65409362767662 +H 1.958764583801193 1.618254016886918 1.67088147423364 0.04710863957682776 -0.03441263562899041 0.09457839337754242 +Li 0.8890679925379519 1.3501094012154145 0.356523142900558 -0.04007135838707443 0.03718283478467768 -0.09125227887537685 +2 +Properties=species:S:1:pos:R:3:forces:R:3 Lattice="2.0 0.0 0.0 0.0 2.0 0.0 0.0 0.0 2.0" energy=-198.27023740453092 +H 1.2841025399739825 1.0560738997199883 0.4076419004771905 -0.6622270842275676 -0.6621719165769338 0.29120726196817465 +Li 0.862900092659493 0.6349124368264709 0.5833565026928345 0.6623678945012492 0.6623120120756423 -0.2882727137301109 +2 +Properties=species:S:1:pos:R:3:forces:R:3 Lattice="2.0 0.0 0.0 0.0 2.0 0.0 0.0 0.0 2.0" energy=-210.62932122368318 +H 1.9990647323447277 0.953558635929386 0.7911340788025603 -0.05594055535694503 -0.057208550417013004 -0.05611647378515233 +Li 1.3290457311009285 0.2507507741141308 0.11711990229764058 0.059840751513080936 0.06045633930206612 0.05994014971060668 +2 +Properties=species:S:1:pos:R:3:forces:R:3 Lattice="2.0 0.0 0.0 0.0 2.0 0.0 0.0 0.0 2.0" energy=-211.11569239713003 +H 0.3293052666603564 0.6745606943807767 1.8083121053115703 0.020429525340863397 -0.02662672320307269 0.005988971418682931 +Li 0.6514315935216397 1.7481427802214673 0.8452546116538187 -0.02444469103096769 0.023429251001441775 0.0019257157545129466 +2 +Properties=species:S:1:pos:R:3:forces:R:3 Lattice="2.0 0.0 0.0 0.0 2.0 0.0 0.0 0.0 2.0" energy=-207.53626571234986 +H 0.8566887767663101 0.514330032046471 0.9820183385163381 -0.23398535394816644 0.194856087268374 0.07061133642254105 +Li 0.10276067711712411 1.0092298537393227 1.167843352890332 0.233035325671203 -0.19619665016142407 -0.07817745982712898 +2 +Properties=species:S:1:pos:R:3:forces:R:3 Lattice="2.0 0.0 0.0 0.0 2.0 0.0 0.0 0.0 2.0" energy=-206.9335273563688 +H 0.9365375039915607 1.070921326767058 1.066203800139247 0.2107657636447698 -0.10299811558112562 0.2770384753981818 +Li 1.3880639850834262 0.836584073664475 1.791586511019425 -0.21262545187531945 0.11296453822579205 -0.27676809894460996 diff --git a/test/data/LiH_crystal_training.xyz b/test/data/LiH_crystal_training.xyz new file mode 100644 index 000000000..ab2791743 --- /dev/null +++ b/test/data/LiH_crystal_training.xyz @@ -0,0 +1,40 @@ +2 +Properties=species:S:1:pos:R:3:forces:R:3 Lattice="2.0 0.0 0.0 0.0 2.0 0.0 0.0 0.0 2.0" energy=-210.34143869968557 +H 1.7469600972641757 1.3105406236661603 1.8078885055514737 -0.05515010202844178 -0.016138638841873887 -0.07494808520830193 +Li 1.3428892907885033 0.34725373984758257 1.1993919048712993 0.05608062426846 0.019533274244299792 0.07536758225127621 +2 +Properties=species:S:1:pos:R:3:forces:R:3 Lattice="2.0 0.0 0.0 0.0 2.0 0.0 0.0 0.0 2.0" energy=-208.1561426455114 +H 0.5421120831942108 1.123751314386705 1.850047836343 0.28176960731737 -0.04513199169175898 0.25950818032265166 +Li 1.228635938071446 1.032454202454703 0.3783686165462776 -0.28071508108709975 0.045080998695914434 -0.26245776947660476 +2 +Properties=species:S:1:pos:R:3:forces:R:3 Lattice="2.0 0.0 0.0 0.0 2.0 0.0 0.0 0.0 2.0" energy=-210.09645903521653 +H 0.027158196041101146 1.6731862097711294 1.6134366515700422 0.10292634370603265 0.1275333232368676 0.04737263331395791 +Li 0.7941335004585335 0.3699666635090202 1.8523215817301748 -0.10497046544138594 -0.12381993344976971 -0.04826856097448054 +2 +Properties=species:S:1:pos:R:3:forces:R:3 Lattice="2.0 0.0 0.0 0.0 2.0 0.0 0.0 0.0 2.0" energy=-205.90144564737625 +H 1.8388327791612642 0.9993481253471506 0.4385729436404704 -0.22101994991201718 -0.16858417104422457 0.34610561623965597 +Li 1.4674006284615426 0.6984981715831435 1.1253342276681788 0.22287617184884612 0.17892796089444607 -0.34832478273018663 +2 +Properties=species:S:1:pos:R:3:forces:R:3 Lattice="2.0 0.0 0.0 0.0 2.0 0.0 0.0 0.0 2.0" energy=-210.78178622831706 +H 1.0168310046613556 0.8612585238423034 0.0009601399098819741 -0.04918730902954849 -0.0494842923887146 -0.002704796073556892 +Li 0.44419313332873056 0.2832776533608796 0.9973760638707188 0.05166094970816282 0.05193243921592955 -0.0008919259574161065 +2 +Properties=species:S:1:pos:R:3:forces:R:3 Lattice="2.0 0.0 0.0 0.0 2.0 0.0 0.0 0.0 2.0" energy=-204.93288634154862 +H 1.6369978211256002 0.3135443027001339 0.17708217034721807 0.286981721245257 0.5591418497105554 -0.3393644337261212 +Li 1.9117219165709483 0.8863594602481111 1.8673689857714335 -0.2855514427745901 -0.5622380627972693 0.33926335540009644 +2 +Properties=species:S:1:pos:R:3:forces:R:3 Lattice="2.0 0.0 0.0 0.0 2.0 0.0 0.0 0.0 2.0" energy=-210.73524383643314 +H 1.7700764233531123 0.7893408914492881 1.088955393401397 -0.048281281265649345 0.05406751470004684 0.0033540969189816416 +Li 1.2528169082220253 1.403038973917562 0.09006815480692909 0.05071086876279346 -0.05608528360044185 0.00022152166059696832 +2 +Properties=species:S:1:pos:R:3:forces:R:3 Lattice="2.0 0.0 0.0 0.0 2.0 0.0 0.0 0.0 2.0" energy=-210.7055029075143 +H 1.6106685007271484 1.0136065810407606 1.3399205019703946 0.07442600324468429 -0.0551332730156715 0.06893938171309927 +Li 0.22790861256500783 0.45451125078373633 0.13344906283841884 -0.07884652763782285 0.05491810155697552 -0.07288355009254194 +2 +Properties=species:S:1:pos:R:3:forces:R:3 Lattice="2.0 0.0 0.0 0.0 2.0 0.0 0.0 0.0 2.0" energy=-211.46640045705647 +H 1.8871080860691987 0.8752458196713906 0.6103866813399352 0.027083053922416 0.0061748501139223255 0.0009190170012606802 +Li 0.7841608913523739 0.008999028007995014 1.3435906437324703 -0.022815401313884376 0.0020283655074712637 -0.009497866107779085 +2 +Properties=species:S:1:pos:R:3:forces:R:3 Lattice="2.0 0.0 0.0 0.0 2.0 0.0 0.0 0.0 2.0" energy=-208.31080335905835 +H 0.751218546267328 1.8293859346784571 0.8411942880456913 0.10733713518674659 -0.1813798156351406 -0.18119888599741052 +Li 1.1181411223414048 1.2045494112405455 0.21726908998803207 -0.11639787090312986 0.18193530292844318 0.18178408033215054 diff --git a/test/test_fchl_acsf_forces.py b/test/test_fchl_acsf_forces.py index 5c870a529..23feeda33 100644 --- a/test/test_fchl_acsf_forces.py +++ b/test/test_fchl_acsf_forces.py @@ -3,7 +3,6 @@ from time import time -import ast import qml from qml.math import cho_solve @@ -36,6 +35,7 @@ np.set_printoptions(linewidth=999, edgeitems=10, suppress=True) import csv +from ase.io import extxyz TRAINING = 7 TEST = 5 @@ -59,6 +59,8 @@ np.random.seed(666) +def MAE(model_output, reference_vals): + return np.mean(np.abs(model_output-reference_vals)) def get_reps(df): @@ -99,6 +101,44 @@ def get_reps(df): return x, f, e, np.array(disp_x), q +def get_reps_pbc(df_pbc): + + x = [] + f_list = [] + e = [] + disp_x = [] + q = [] + + for i in df_pbc: + + coordinates = i.get_positions() + nuclear_charges = i.get_atomic_numbers() + atomtypes = [1, 3] + + force = np.array(i.get_forces()) + energy = float(i.get_total_energy()) + cell = i.cell[:] + + (x1, dx1) = generate_fchl_acsf(nuclear_charges, coordinates, gradients=True, \ + rcut=1.95, elements = atomtypes, cell = cell) + + x.append(x1) + f_list.append(force) + e.append(energy) + + disp_x.append(dx1) + q.append(nuclear_charges) + + e = np.array(e) + # e -= np.mean(e)# - 10 # + + f = np.stack(f_list) + f *= -1 + x = np.array(x) + + return x, f, e, np.array(disp_x), q + + def test_fchl_acsf_operator(): print("Representations ...") @@ -165,6 +205,60 @@ def test_fchl_acsf_operator(): (np.mean(np.abs(Fv.flatten() - fYv.flatten())), slope, intercept, r_value )) +def test_fchl_acsf_operator_pbc(): + + SIGMA = 20.0 + TRAINING = 10 + TEST = 8 + + + print("Representations (pbc) ...") + + xyz_pbc_training = open(TEST_DIR+"/data/LiH_crystal_training.xyz", 'r') + DF_TRAIN_PBC = list(extxyz.read_extxyz(xyz_pbc_training, index=slice(None, None, 1))) + xyz_pbc_training.close() + + xyz_pbc_test = open(TEST_DIR+"/data/LiH_crystal_test.xyz", 'r') + DF_TEST_PBC = list(extxyz.read_extxyz(xyz_pbc_test, index=slice(None, None, 1))) + xyz_pbc_test.close() + + X, F, E, dX, Q = get_reps_pbc(DF_TRAIN_PBC) + Xs, Fs, Es, dXs, Qs = get_reps_pbc(DF_TEST_PBC) + + F = np.concatenate(F) + Fs = np.concatenate(Fs) + + Y = np.concatenate((E, F.flatten())) + + print("Kernels (pbc) ...") + Kte = get_atomic_local_kernel(X, X, Q, Q, SIGMA) + Kse = get_atomic_local_kernel(X, Xs, Q, Qs, SIGMA) + + Kt = get_atomic_local_gradient_kernel(X, X, dX, Q, Q, SIGMA) + Ks = get_atomic_local_gradient_kernel(X, Xs, dXs, Q, Qs, SIGMA) + + C = np.concatenate((Kte, Kt)) + + print("Alphas operator (pbc) ...") + alpha = svd_solve(C, Y, rcond=1e-9) + # alpha = qrlq_solve(C, Y) + + eYt = np.dot(Kte, alpha) + eYs = np.dot(Kse, alpha) + + fYt = np.dot(Kt, alpha) + fYs = np.dot(Ks, alpha) + + + print("===============================================================================================") + print("==== OPERATOR, FORCE + ENERGY (PBC) ==========================================================") + print("===============================================================================================") + + print("TEST ENERGY MAE = %10.6f MAE (expected) = %10.6f " % (np.mean(np.abs(Es - eYs)), 2.363580)) + print("TEST FORCE MAE = %10.6f MAE (expected) = %10.6f " % (np.mean(np.abs(Fs.flatten() - fYs.flatten())), 0.981332)) + + + def test_fchl_acsf_gaussian_process(): print("Representations ...") From 68ef9ea3d7292be5713004114df92ef5f9ef12a4 Mon Sep 17 00:00:00 2001 From: kvkarandashev Date: Thu, 4 Mar 2021 19:32:30 +0100 Subject: [PATCH 3/8] Updated files' headers. --- qml/kernels/fgradient_kernels.f90 | 2 +- qml/representations/facsf.f90 | 2 +- qml/representations/representations.py | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/qml/kernels/fgradient_kernels.f90 b/qml/kernels/fgradient_kernels.f90 index 5bdb45f29..55f9d13de 100644 --- a/qml/kernels/fgradient_kernels.f90 +++ b/qml/kernels/fgradient_kernels.f90 @@ -1,6 +1,6 @@ ! MIT License ! -! Copyright (c) 2018-2019 Anders Steen Christensen +! Copyright (c) 2018-2021 Anders Steen Christensen, Konstantin Karandashev ! ! Permission is hereby granted, free of charge, to any person obtaining a copy ! of this software and associated documentation files (the "Software"), to deal diff --git a/qml/representations/facsf.f90 b/qml/representations/facsf.f90 index 25e587a98..230f787ea 100644 --- a/qml/representations/facsf.f90 +++ b/qml/representations/facsf.f90 @@ -1,6 +1,6 @@ ! MIT License ! -! Copyright (c) 2018 Lars Andersen Bratholm +! Copyright (c) 2018-2021 Lars Andersen Bratholm, Konstantin Karandashev ! ! Permission is hereby granted, free of charge, to any person obtaining a copy ! of this software and associated documentation files (the "Software"), to deal diff --git a/qml/representations/representations.py b/qml/representations/representations.py index a8f2b5386..f3f9be9ed 100644 --- a/qml/representations/representations.py +++ b/qml/representations/representations.py @@ -1,6 +1,6 @@ # MIT License # -# Copyright (c) 2017-2018 Anders Steen Christensen, Lars Andersen Bratholm, Bing Huang +# Copyright (c) 2017-2021 Anders Steen Christensen, Lars Andersen Bratholm, Bing Huang, Konstantin Karandashev # # Permission is hereby granted, free of charge, to any person obtaining a copy # of this software and associated documentation files (the "Software"), to deal From 55afd825aec0c9620ef5b99e828fe3a496b89191 Mon Sep 17 00:00:00 2001 From: kvkarandashev Date: Tue, 9 Mar 2021 18:02:53 +0100 Subject: [PATCH 4/8] Second draft for a test for PBC implementation for FCHL19. --- test/test_fchl_acsf_pbc.py | 120 +++++++++++++++++++++++++++++++++++++ 1 file changed, 120 insertions(+) create mode 100644 test/test_fchl_acsf_pbc.py diff --git a/test/test_fchl_acsf_pbc.py b/test/test_fchl_acsf_pbc.py new file mode 100644 index 000000000..52c5ff102 --- /dev/null +++ b/test/test_fchl_acsf_pbc.py @@ -0,0 +1,120 @@ +# MIT License +# +# Copyright (c) 2021 Konstantin Karandashev +# +# Permission is hereby granted, free of charge, to any person obtaining a copy +# of this software and associated documentation files (the "Software"), to deal +# in the Software without restriction, including without limitation the rights +# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +# copies of the Software, and to permit persons to whom the Software is +# furnished to do so, subject to the following conditions: +# +# The above copyright notice and this permission notice shall be included in all +# copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +# SOFTWARE. + +""" +Runs the same calculation using the 'cell' keyword for generate_fchl_acsf and by straightforward +'cell cloning' done inside the test script. +""" +import os +from copy import deepcopy +import numpy as np +np.set_printoptions(linewidth=666, edgeitems=10) +from qml import Compound +from qml.representations import generate_fchl_acsf +import random + +REP_PARAMS = dict() +REP_PARAMS["elements"] = [1, 6, 7, 8, 16] +REP_PARAMS["rcut"]=5.0 +random.seed(1) +cell_added_cutoff=0.1 + +# For given molecular coordinates generate a cell just large enough to contain the molecule. +def suitable_cell(coords): + max_coords=None + min_coords=None + for atom_coords in coords: + if max_coords is None: + max_coords=deepcopy(atom_coords) + min_coords=deepcopy(atom_coords) + else: + max_coords=np.maximum(max_coords, atom_coords) + min_coords=np.minimum(min_coords, atom_coords) + return np.diag((max_coords-min_coords)*(1.0+cell_added_cutoff)) + +def generate_fchl_acsf_brute_pbc(nuclear_charges, coordinates, cell, gradients=False): + num_atoms=len(nuclear_charges) + all_coords=deepcopy(coordinates) + all_charges=deepcopy(nuclear_charges) + nExtend = (np.floor(REP_PARAMS["rcut"]/np.linalg.norm(cell,2,axis = 0)) + 1).astype(int) + print("Checked nExtend:", nExtend, ", gradient calculation:", gradients) + for i in range(-nExtend[0],nExtend[0] + 1): + for j in range(-nExtend[1],nExtend[1] + 1): + for k in range(-nExtend[2],nExtend[2] + 1): + if not (i == 0 and j == 0 and k == 0): + all_coords=np.append(all_coords, coordinates + i*cell[0,:] + j*cell[1,:] + k*cell[2,:], axis=0) + all_charges=np.append(all_charges, nuclear_charges) + if gradients: + rep, drep=generate_fchl_acsf(all_charges, all_coords, gradients=True, **REP_PARAMS) + else: + rep=generate_fchl_acsf(all_charges, all_coords, gradients=False, **REP_PARAMS) + rep=rep[:num_atoms,:] + if gradients: + drep=drep[:num_atoms,:,:num_atoms,:] + return rep, drep + else: + return rep + +def ragged_array_close(arr1, arr2, error_msg): + for el1, el2 in zip(arr1, arr2): + assert np.allclose(el1, el2), error_msg + +def test_fchl_acsf_pbc(): + + qm7_dir = os.path.dirname(os.path.realpath(__file__))+"/qm7" + os.chdir(qm7_dir) + all_xyzs=os.listdir() + test_xyzs=random.sample(all_xyzs, 3) + + reps_no_grad1=[] + reps_no_grad2=[] + + reps_wgrad1=[] + reps_wgrad2=[] + dreps1=[] + dreps2=[] + + for xyz in test_xyzs: + print("Tested xyz:", xyz) + mol = Compound(xyz=xyz) + cell=suitable_cell(mol.coordinates) + reps_no_grad1.append(generate_fchl_acsf_brute_pbc(mol.nuclear_charges, mol.coordinates, cell, gradients=False)) + reps_no_grad2.append(generate_fchl_acsf(mol.nuclear_charges, mol.coordinates, cell=cell, gradients=False, **REP_PARAMS)) + + rep_wgrad1, drep1=generate_fchl_acsf_brute_pbc(mol.nuclear_charges, mol.coordinates, cell, gradients=True) + rep_wgrad2, drep2=generate_fchl_acsf(mol.nuclear_charges, mol.coordinates, cell=cell, gradients=True, **REP_PARAMS) + + reps_wgrad1.append(rep_wgrad1) + reps_wgrad2.append(rep_wgrad2) + + dreps1.append(drep1) + dreps2.append(drep2) + + ragged_array_close(reps_no_grad1, reps_no_grad2, "Error in PBC implementation for generate_fchl_acsf without gradients.") + ragged_array_close(reps_wgrad1, reps_wgrad2, "Error in PBC implementation for generate_fchl_acsf with gradients (representation).") + ragged_array_close(dreps1, dreps2, "Error in PBC implementation for generate_fchl_acsf with gradients (gradient of representation).") + print("Passed") + +if __name__ == "__main__": + + test_fchl_acsf_pbc() + From 9faadbccc27ed85b5a160d4f274f0e03a3ea8ce2 Mon Sep 17 00:00:00 2001 From: kvkarandashev Date: Mon, 1 Apr 2024 16:14:22 +0200 Subject: [PATCH 5/8] Minor FCHL revision. --- qml/representations/facsf.f90 | 328 ++++++++++++------------- qml/representations/representations.py | 2 +- test/test_fchl_acsf_pbc.py | 5 +- 3 files changed, 166 insertions(+), 169 deletions(-) diff --git a/qml/representations/facsf.f90 b/qml/representations/facsf.f90 index 230f787ea..a918b1793 100644 --- a/qml/representations/facsf.f90 +++ b/qml/representations/facsf.f90 @@ -1,24 +1,24 @@ -! MIT License -! -! Copyright (c) 2018-2021 Lars Andersen Bratholm, Konstantin Karandashev -! -! Permission is hereby granted, free of charge, to any person obtaining a copy -! of this software and associated documentation files (the "Software"), to deal -! in the Software without restriction, including without limitation the rights -! to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -! copies of the Software, and to permit persons to whom the Software is -! furnished to do so, subject to the following conditions: -! -! The above copyright notice and this permission notice shall be included in all -! copies or substantial portions of the Software. -! -! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -! IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -! AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -! SOFTWARE. + + + + + + + + + + + + + + + + + + + + + module acsf_utils @@ -434,13 +434,13 @@ subroutine fgenerate_acsf_and_gradients(coordinates, nuclear_charges, elements, do k = 1, 3 ! The gradients wrt coordinates part = radial_part * (coordinates(i,k) - coordinates(j,k)) - grad_subset(i, (n-1)*nbasis2 + 1:n*nbasis2, i, k) = & + grad_subset(i, (n-1)*nbasis2 + 1:n*nbasis2, i, k) = & grad_subset(i, (n-1)*nbasis2 + 1:n*nbasis2, i, k) + part - grad_subset(i, (n-1)*nbasis2 + 1:n*nbasis2, j, k) = & + grad_subset(i, (n-1)*nbasis2 + 1:n*nbasis2, j, k) = & grad_subset(i, (n-1)*nbasis2 + 1:n*nbasis2, j, k) - part - grad_subset(j, (m-1)*nbasis2 + 1:m*nbasis2, j, k) = & + grad_subset(j, (m-1)*nbasis2 + 1:m*nbasis2, j, k) = & grad_subset(j, (m-1)*nbasis2 + 1:m*nbasis2, j, k) - part - grad_subset(j, (m-1)*nbasis2 + 1:m*nbasis2, i, k) = & + grad_subset(j, (m-1)*nbasis2 + 1:m*nbasis2, i, k) = & grad_subset(j, (m-1)*nbasis2 + 1:m*nbasis2, i, k) + part enddo endif @@ -593,7 +593,7 @@ subroutine fgenerate_acsf_and_gradients(coordinates, nuclear_charges, elements, atom_grad(z:z + nabasis - 1, k, t) = atom_grad(z:z + nabasis - 1, k, t) + & & d_angular * d_angular_d_k(t) * radial(l) + & & angular * d_radial(l) * d_radial_d_k(t) - & - & angular * radial(l) * rdecay(i,j) * d_ikdecay(t) + & angular * radial(l) * rdecay(i,j) * d_ikdecay(t) enddo enddo enddo @@ -659,14 +659,16 @@ subroutine fgenerate_fchl_acsf(coordinates, nuclear_charges, elements, & double precision, intent(in) :: three_body_weight double precision, intent(out), dimension(natoms, rep_size) :: rep -! Introduced to make OpenMP parallelization easier. - double precision, dimension(:, :, :), allocatable:: add_rep - integer:: i, j, k, l, n, m, o, p, q, s, z, nelements, nbasis2, nbasis3, nabasis + integer:: two_body_rep_size + integer :: i, j, k, l, n, m, o, p, q, s, z, nelements, nbasis2, nbasis3, nabasis, j_id, j_id_max integer, allocatable, dimension(:) :: element_types double precision :: rij, rik, angle, cos_1, cos_2, cos_3, invcut + ! double precision :: angle_1, angle_2, angle_3 double precision, allocatable, dimension(:) :: radial, angular, a, b, c double precision, allocatable, dimension(:, :) :: distance_matrix, rdecay + double precision, allocatable, dimension(:, :) :: saved_radial + integer, allocatable, dimension(:):: saved_j double precision :: mu, sigma, ksi3 @@ -687,7 +689,7 @@ subroutine fgenerate_fchl_acsf(coordinates, nuclear_charges, elements, & allocate(element_types(natoms_tot)) ! Store element index of every atom -!$OMP PARALLEL DO SCHEDULE(dynamic) + !$OMP PARALLEL DO SCHEDULE(dynamic) do i = 1, natoms_tot do j = 1, nelements if (nuclear_charges(modulo(i-1, natoms)+1) .eq. elements(j)) then @@ -696,7 +698,7 @@ subroutine fgenerate_fchl_acsf(coordinates, nuclear_charges, elements, & endif enddo enddo -!$OMP END PARALLEL DO + !$OMP END PARALLEL DO ! Get distance matrix @@ -705,7 +707,7 @@ subroutine fgenerate_fchl_acsf(coordinates, nuclear_charges, elements, & distance_matrix = 0.0d0 -!$OMP PARALLEL DO PRIVATE(rij) SCHEDULE(dynamic) + !$OMP PARALLEL DO PRIVATE(rij) SCHEDULE(dynamic) do i = 1, natoms_tot do j = i+1, natoms_tot rij = norm2(coordinates(j,:) - coordinates(i,:)) @@ -713,7 +715,7 @@ subroutine fgenerate_fchl_acsf(coordinates, nuclear_charges, elements, & distance_matrix(j, i) = rij enddo enddo -!$OMP END PARALLEL DO + !$OMP END PARALLEL DO ! number of basis functions in the two body term nbasis2 = size(Rs2) @@ -725,49 +727,46 @@ subroutine fgenerate_fchl_acsf(coordinates, nuclear_charges, elements, & rdecay = decay(distance_matrix, invcut, natoms_tot) ! Allocate temporary - allocate(radial(nbasis2), add_rep(nbasis2, natoms, natoms)) - - radial = 0.0d0 - add_rep=0.0d0 -!$OMP PARALLEL DO PRIVATE(n,m,rij,radial,mu,sigma) SCHEDULE(dynamic) + allocate(radial(nbasis2)) + allocate(saved_radial(nbasis2, natoms_tot), saved_j(natoms_tot)) + rep=0.0d0 + !$OMP PARALLEL DO PRIVATE(n,m,rij,radial,mu,sigma,saved_radial, saved_j, j_id_max) SCHEDULE(dynamic) do i = 1, natoms ! index of the element of atom i m = element_types(i) + j_id_max=0 do j = i + 1, natoms_tot - ! index of the element of atom j - n = element_types(j) ! distance between atoms i and j rij = distance_matrix(i,j) - if (rij <= rcut) then - - ! two body term of the representation - mu = log(rij / sqrt(1.0d0 + eta2 / rij**2)) - sigma = sqrt(log(1.0d0 + eta2 / rij**2)) - radial(:) = 0.0d0 + if (rij > rcut) cycle + j_id_max=j_id_max+1 + saved_j(j_id_max)=j + ! index of the element of atom j + n = element_types(j) + ! two body term of the representation + mu = log(rij / sqrt(1.0d0 + eta2 / rij**2)) + sigma = sqrt(log(1.0d0 + eta2 / rij**2)) + radial(:) = 0.0d0 - do k = 1, nbasis2 - radial(k) = 1.0d0/(sigma* sqrt(2.0d0*pi) * Rs2(k)) * rdecay(i,j) & + do k = 1, nbasis2 + radial(k) = 1.0d0/(sigma* sqrt(2.0d0*pi) * Rs2(k)) * rdecay(i,j) & & * exp( - (log(Rs2(k)) - mu)**2 / (2.0d0 * sigma**2) ) / rij**two_body_decay - enddo - - rep(i, (n-1)*nbasis2 + 1:n*nbasis2) = rep(i, (n-1)*nbasis2 + 1:n*nbasis2) + radial -! rep(j, (m-1)*nbasis2 + 1:m*nbasis2) = rep(j, (m-1)*nbasis2 + 1:m*nbasis2) + radial - if (j<=natoms) add_rep(:, i, j)=radial - endif + enddo + saved_radial(:, j_id_max)=radial(:) enddo - enddo -!$OMP END PARALLEL DO - deallocate(radial) -!$OMP PARALLEL DO PRIVATE(m) SCHEDULE(dynamic) - do j = 2, natoms - do i = 1, j-1 - m = element_types(i) - rep(j, (m-1)*nbasis2 + 1:m*nbasis2) = rep(j, (m-1)*nbasis2 + 1:m*nbasis2) + add_rep(:, i, j) + !$OMP CRITICAL + do j_id = 1, j_id_max + j=saved_j(j_id) + n=element_types(j) + rep(i, (n-1)*nbasis2 + 1:n*nbasis2) = rep(i, (n-1)*nbasis2 + 1:n*nbasis2) + saved_radial(:, j_id) + if (j<=natoms) rep(j, (m-1)*nbasis2 + 1:m*nbasis2) = rep(j, (m-1)*nbasis2 + 1:m*nbasis2) + saved_radial(:, j_id) enddo + !$OMP END CRITICAL enddo -!$OMP END PARALLEL DO + !$OMP END PARALLEL DO - deallocate(add_rep) + deallocate(radial) + deallocate(saved_radial, saved_j) ! number of radial basis functions in the three body term nbasis3 = size(Rs3) @@ -785,9 +784,13 @@ subroutine fgenerate_fchl_acsf(coordinates, nuclear_charges, elements, & allocate(c(3)) allocate(radial(nbasis3)) allocate(angular(nabasis)) + two_body_rep_size=nbasis2*nelements -!$OMP PARALLEL DO PRIVATE(rij, n, rik, m, a, b, c, angle, cos_1, cos_2, cos_3,& -!$OMP radial, ksi3, angular, o, p, q, s, z) SCHEDULE(dynamic) + ! This could probably be done more efficiently if it's a bottleneck + ! Also the order is a bit wobbly compared to the tensorflow implementation + !$OMP PARALLEL DO PRIVATE(rij, n, rik, m, a, b, c, angle, radial, angular, & + !$OMP cos_1, cos_2, cos_3, mu, sigma, o, ksi3, & + !$OMP p, q, s, z, l) SCHEDULE(dynamic) do i = 1, natoms do j = 1, natoms_tot - 1 if (i .eq. j) cycle @@ -816,29 +819,29 @@ subroutine fgenerate_fchl_acsf(coordinates, nuclear_charges, elements, & ! The radial part of the three body terms including decay radial = exp(-eta3*(0.5d0 * (rij+rik) - Rs3)**2) * rdecay(i,j) * rdecay(i,k) - + ksi3 = (1.0d0 + 3.0d0 * cos_1 * cos_2 * cos_3) & & / (distance_matrix(i,k) * distance_matrix(i,j) * distance_matrix(j,k) & & )**three_body_decay * three_body_weight - angular = 0.0d0 + angular = 0.0d0 do l = 1, nabasis/2 o = l*2-1 angular(2*l-1) = angular(2*l-1) + 2*cos(o * angle) & & * exp(-(zeta * o)**2 /2) - + angular(2*l) = angular(2*l) + 2*sin(o * angle) & & * exp(-(zeta * o)**2 /2) enddo - + ! The lowest of the element indices for atoms j and k p = min(n,m) - 1 ! The highest of the element indices for atoms j and k q = max(n,m) - 1 ! calculate the indices that the three body terms should be added to - s = nelements * nbasis2 + nbasis3 * nabasis * (-(p * (p + 1))/2 + q + nelements * p) + 1 + s = two_body_rep_size + nbasis3 * nabasis * (-(p * (p + 1))/2 + q + nelements * p) + 1 do l = 1, nbasis3 ! calculate the indices that the three body terms should be added to @@ -849,7 +852,7 @@ subroutine fgenerate_fchl_acsf(coordinates, nuclear_charges, elements, & enddo enddo enddo -!$OMP END PARALLEL DO + !$OMP END PARALLEL DO deallocate(element_types) deallocate(rdecay) @@ -863,8 +866,8 @@ subroutine fgenerate_fchl_acsf(coordinates, nuclear_charges, elements, & end subroutine fgenerate_fchl_acsf subroutine fgenerate_fchl_acsf_and_gradients(coordinates, nuclear_charges, elements, & - & Rs2, Rs3, Ts, eta2, eta3, zeta, rcut, acut, natoms, natoms_tot, rep_size, & - & two_body_decay, three_body_decay, three_body_weight, rep, grad) + & Rs2, Rs3, Ts, eta2, eta3, zeta, rcut, acut, natoms, natoms_tot, rep_size, & + & two_body_decay, three_body_decay, three_body_weight, rep, grad) use acsf_utils, only: decay, calc_angle, calc_cos_angle @@ -881,7 +884,7 @@ subroutine fgenerate_fchl_acsf_and_gradients(coordinates, nuclear_charges, eleme double precision, intent(in) :: zeta double precision, intent(in) :: rcut double precision, intent(in) :: acut - + double precision, intent(in) :: two_body_decay double precision, intent(in) :: three_body_decay double precision, intent(in) :: three_body_weight @@ -927,18 +930,17 @@ subroutine fgenerate_fchl_acsf_and_gradients(coordinates, nuclear_charges, eleme if (natoms /= size(nuclear_charges, dim=1)) then write(*,*) "ERROR: Atom Centered Symmetry Functions creation" write(*,*) natoms, "coordinates, but", & - & size(nuclear_charges, dim=1), "atom_types!" + & size(nuclear_charges, dim=1), "atom_types!" stop endif - ! Number of unique elements nelements = size(elements) ! Allocate temporary allocate(element_types(natoms_tot)) -! Store element index of every atom -!$OMP PARALLEL DO SCHEDULE(dynamic) + ! Store element index of every atom + !$OMP PARALLEL DO SCHEDULE(dynamic) do i = 1, natoms_tot do j = 1, nelements if (nuclear_charges(modulo(i-1, natoms)+1) .eq. elements(j)) then @@ -947,7 +949,7 @@ subroutine fgenerate_fchl_acsf_and_gradients(coordinates, nuclear_charges, eleme endif enddo enddo -!$OMP END PARALLEL DO + !$OMP END PARALLEL DO @@ -963,7 +965,7 @@ subroutine fgenerate_fchl_acsf_and_gradients(coordinates, nuclear_charges, eleme inv_sq_distance_matrix = 0.0d0 -!$OMP PARALLEL DO PRIVATE(rij,rij2,invrij,invrij2) SCHEDULE(dynamic) + !$OMP PARALLEL DO PRIVATE(rij,rij2,invrij,invrij2) SCHEDULE(dynamic) do i = 1, natoms_tot do j = i+1, natoms_tot rij = norm2(coordinates(j,:) - coordinates(i,:)) @@ -980,7 +982,7 @@ subroutine fgenerate_fchl_acsf_and_gradients(coordinates, nuclear_charges, eleme inv_sq_distance_matrix(j, i) = invrij2 enddo enddo -!$OMP END PARALLEL DO + !$OMP END PARALLEL DO ! Number of two body basis functions @@ -1007,75 +1009,69 @@ subroutine fgenerate_fchl_acsf_and_gradients(coordinates, nuclear_charges, eleme add_grad=0.0d0 log_Rs2(:) = log(Rs2(:)) -!$OMP PARALLEL DO PRIVATE(m, n, rij, invrij, mu, sigma, exp_s2, exp_ln,& -!$OMP scaling, radial_base, radial, dx, part, dscal, ddecay) SCHEDULE(dynamic) + !$OMP PARALLEL DO PRIVATE(m, n, rij, invrij, mu, sigma, exp_s2, exp_ln,& + !$OMP scaling, radial_base, radial, dx, part, dscal, ddecay) SCHEDULE(dynamic) do i = 1, natoms ! The element index of atom i - m = element_types(i) + m = element_types(i) ! moved inside loop to enable COLLAPSE do j = i + 1, natoms_tot + rij = distance_matrix(i,j) + if (rij > rcut) continue ! The element index of atom j n = element_types(j) ! Distance between atoms i and j - rij = distance_matrix(i,j) - if (rij <= rcut) then - invrij = inv_distance_matrix(i,j) - - mu = log(rij / sqrt(1.0d0 + eta2 * inv_sq_distance_matrix(i, j))) - sigma = sqrt(log(1.0d0 + eta2 * inv_sq_distance_matrix(i, j))) - exp_s2 = exp(sigma**2) - exp_ln = exp(-(log_Rs2(:) - mu)**2 / sigma**2 * 0.5d0) * sqrt(2.0d0) - - scaling = 1.0d0 / rij**two_body_decay + invrij = inv_distance_matrix(i,j) + mu = log(rij / sqrt(1.0d0 + eta2 * inv_sq_distance_matrix(i, j))) + sigma = sqrt(log(1.0d0 + eta2 * inv_sq_distance_matrix(i, j))) + exp_s2 = exp(sigma**2) + exp_ln = exp(-(log_Rs2(:) - mu)**2 / sigma**2 * 0.5d0) * sqrt(2.0d0) + scaling = 1.0d0 / rij**two_body_decay - radial_base(:) = 1.0d0/(sigma* sqrt(2.0d0*pi) * Rs2(:)) * exp(-(log_Rs2(:) - mu)**2 / (2.0d0 * sigma**2)) - radial(:) = radial_base(:) * scaling * rdecay(i,j) - - rep(i, (n-1)*nbasis2 + 1:n*nbasis2) = rep(i, (n-1)*nbasis2 + 1:n*nbasis2) + radial -! rep(j, (m-1)*nbasis2 + 1:m*nbasis2) = rep(j, (m-1)*nbasis2 + 1:m*nbasis2) + radial - if (j<=natoms) add_rep(:, i, j)=radial + radial_base(:) = 1.0d0/(sigma* sqrt(2.0d0*pi) * Rs2(:)) * exp(-(log_Rs2(:) - mu)**2 / (2.0d0 * sigma**2)) - do k = 1, 3 + radial(:) = radial_base(:) * scaling * rdecay(i,j) - dx = -(coordinates(i,k) - coordinates(j,k)) - - part(:) = ((log_Rs2(:) - mu) * (-dx *(rij**2 * exp_s2 + eta2) / (rij * sqrt(exp_s2))**3) & - &* sqrt(exp_s2) / (sigma**2 * rij) + (log_Rs2(:) - mu) ** 2 * eta2 * dx / & - &(sigma**4 * rij**4 * exp_s2)) * exp_ln / (Rs2(:) * sigma * sqrt(pi) * 2) & - &- exp_ln * eta2 * dx / (Rs2(:) * sigma**3 *sqrt(pi) * rij**4 * exp_s2 * 2.0d0) + rep(i, (n-1)*nbasis2 + 1:n*nbasis2) = rep(i, (n-1)*nbasis2 + 1:n*nbasis2) + radial + if (j<=natoms) add_rep(:, i, j)=radial + do k = 1, 3 + dx = -(coordinates(i,k) - coordinates(j,k)) + part(:) = ((log_Rs2(:) - mu) * (-dx *(rij**2 * exp_s2 + eta2) / (rij * sqrt(exp_s2))**3) & + &* sqrt(exp_s2) / (sigma**2 * rij) + (log_Rs2(:) - mu) ** 2 * eta2 * dx / & + &(sigma**4 * rij**4 * exp_s2)) * exp_ln / (Rs2(:) * sigma * sqrt(pi) * 2) & + &- exp_ln * eta2 * dx / (Rs2(:) * sigma**3 *sqrt(pi) * rij**4 * exp_s2 * 2.0d0) - dscal = two_body_decay * dx / rij**(two_body_decay+2.0d0) - ddecay = dx * 0.5d0 * pi * sin(pi*rij * invcut) * invcut * invrij + dscal = two_body_decay * dx / rij**(two_body_decay+2.0d0) + ddecay = dx * 0.5d0 * pi * sin(pi*rij * invcut) * invcut * invrij - part(:) = part(:) * scaling * rdecay(i,j) + radial_base(:) * dscal * rdecay(i,j) & - & + radial_base(:) * scaling * ddecay + part(:) = part(:) * scaling * rdecay(i,j) + radial_base(:) * dscal * rdecay(i,j) & + & + radial_base(:) * scaling * ddecay - ! The gradients wrt coordinates - grad(i, (n-1)*nbasis2 + 1:n*nbasis2, i, k) = grad(i, (n-1)*nbasis2 + 1:n*nbasis2, i, k) + part - if (j<=natoms) then - grad(i, (n-1)*nbasis2 + 1:n*nbasis2, j, k) = grad(i, (n-1)*nbasis2 + 1:n*nbasis2, j, k) - part - grad(j, (m-1)*nbasis2 + 1:m*nbasis2, i, k) = grad(j, (m-1)*nbasis2 + 1:m*nbasis2, i, k) + part -! grad(j, (m-1)*nbasis2 + 1:m*nbasis2, j, k) = grad(j, (m-1)*nbasis2 + 1:m*nbasis2, j, k) - part - add_grad(:, k, i, j)=part - endif - enddo - endif + ! The gradients wrt coordinates + grad(i, (n-1)*nbasis2 + 1:n*nbasis2, i, k) = grad(i, (n-1)*nbasis2 + 1:n*nbasis2, i, k) + part + if (j<=natoms) then + grad(i, (n-1)*nbasis2 + 1:n*nbasis2, j, k) = grad(i, (n-1)*nbasis2 + 1:n*nbasis2, j, k) - part + grad(j, (m-1)*nbasis2 + 1:m*nbasis2, i, k) = grad(j, (m-1)*nbasis2 + 1:m*nbasis2, i, k) + part + add_grad(:, k, i, j)=part + endif + enddo enddo enddo -!$OMP END PARALLEL DO + !$OMP END PARALLEL DO -!$OMP PARALLEL DO PRIVATE(m) SCHEDULE(dynamic) + !$OMP PARALLEL DO PRIVATE(m) SCHEDULE(dynamic) do j = 2, natoms do i = 1, j-1 m = element_types(i) + if (distance_matrix(i,j)>rcut) cycle rep(j, (m-1)*nbasis2 + 1:m*nbasis2) = rep(j, (m-1)*nbasis2 + 1:m*nbasis2) + add_rep(:, i, j) do k=1, 3 grad(j, (m-1)*nbasis2 + 1:m*nbasis2, j, k) = grad(j, (m-1)*nbasis2 + 1:m*nbasis2, j, k) - add_grad(:, k, i, j) - enddo + enddo enddo enddo -!$OMP END PARALLEL DO + !$OMP END PARALLEL DO deallocate(add_rep, add_grad) @@ -1140,13 +1136,13 @@ subroutine fgenerate_fchl_acsf_and_gradients(coordinates, nuclear_charges, eleme allocate(d_atm_extra_j(3)) allocate(d_atm_extra_k(3)) -!$OMP PARALLEL DO PRIVATE(atom_rep, atom_grad, rij, n, rij2, invrij, invrij2,& -!$OMP rik, m, rik2, invrik, invrjk, invrik2, a, b, c, angle, cos_i, cos_k,& -!$OMP cos_j, radial_base, radial, p, q, dot, angular, d_angular, d_angular_d_j,& -!$OMP d_angular_d_k, d_angular_d_i, d_radial, d_radial_d_j, d_radial_d_k,& -!$OMP d_radial_d_i, d_ijdecay, d_ikdecay, invr_atm, atm, atm_i, atm_j, atm_k, vi,& -!$OMP vj, vk, d_atm_ii, d_atm_ij, d_atm_ik, d_atm_ji, d_atm_jj, d_atm_jk, d_atm_ki,& -!$OMP d_atm_kj, d_atm_kk, d_atm_extra_i, d_atm_extra_j, d_atm_extra_k, s, z) SCHEDULE(dynamic) + !$OMP PARALLEL DO PRIVATE(atom_rep, atom_grad, rij, n, rij2, invrij, invrij2,& + !$OMP rik, m, rik2, invrik, invrjk, invrik2, a, b, c, angle, cos_i, cos_k,& + !$OMP cos_j, radial_base, radial, p, q, dot, angular, d_angular, d_angular_d_j,& + !$OMP d_angular_d_k, d_angular_d_i, d_radial, d_radial_d_j, d_radial_d_k,& + !$OMP d_radial_d_i, d_ijdecay, d_ikdecay, invr_atm, atm, atm_i, atm_j, atm_k, vi,& + !$OMP vj, vk, d_atm_ii, d_atm_ij, d_atm_ik, d_atm_ji, d_atm_jj, d_atm_jk, d_atm_ki,& + !$OMP d_atm_kj, d_atm_kk, d_atm_extra_i, d_atm_extra_j, d_atm_extra_k, s, z) SCHEDULE(dynamic) do i = 1, natoms atom_rep = 0.0d0 atom_grad = 0.0d0 @@ -1187,7 +1183,7 @@ subroutine fgenerate_fchl_acsf_and_gradients(coordinates, nuclear_charges, eleme cos_i = calc_cos_angle(a,b,c) cos_k = calc_cos_angle(a,c,b) cos_j = calc_cos_angle(b,a,c) - + ! part of the radial part of the 3body terms radial_base(:) = exp(-eta3*(0.5d0 * (rij+rik) - Rs3(:))**2) radial(:) = radial_base(:) ! * scaling @@ -1197,7 +1193,7 @@ subroutine fgenerate_fchl_acsf_and_gradients(coordinates, nuclear_charges, eleme q = max(n,m) - 1 ! Dot product between the vectors connecting atom i,j and i,k dot = dot_product(a-b,c-b) - + angular(1) = exp(-(zeta**2)*0.5d0) * 2 * cos(angle) angular(2) = exp(-(zeta**2)*0.5d0) * 2 * sin(angle) @@ -1225,7 +1221,7 @@ subroutine fgenerate_fchl_acsf_and_gradients(coordinates, nuclear_charges, eleme d_ijdecay = - pi * (b - a) * sin(pi * rij * invcut) * 0.5d0 * invrij * invcut ! Part of the derivative of the i,k decay functions wrt coordinates (dim(3)) d_ikdecay = - pi * (b - c) * sin(pi * rik * invcut) * 0.5d0 * invrik * invcut - + invr_atm = (invrij * invrjk *invrik)**three_body_decay ! Axilrod-Teller-Muto term @@ -1234,15 +1230,15 @@ subroutine fgenerate_fchl_acsf_and_gradients(coordinates, nuclear_charges, eleme atm_i = (3.0d0 * cos_j * cos_k) * invr_atm * invrij * invrik atm_j = (3.0d0 * cos_k * cos_i) * invr_atm * invrij * invrjk atm_k = (3.0d0 * cos_i * cos_j) * invr_atm * invrjk * invrik - + vi = dot_product(a-b,c-b) vj = dot_product(c-a,b-a) vk = dot_product(b-c,a-c) - + d_atm_ii(:) = 2 * b - a - c - vi * ((b-a)*invrij**2 + (b-c)*invrik**2) d_atm_ij(:) = c - a - vj * (b-a)*invrij**2 d_atm_ik(:) = a - c - vk * (b-c)*invrik**2 - + d_atm_ji(:) = c - b - vi * (a-b)*invrij**2 d_atm_jj(:) = 2 * a - b - c - vj * ((a-b)*invrij**2 + (a-c)*invrjk**2) d_atm_jk(:) = b - c - vk * (a-c)*invrjk**2 @@ -1268,7 +1264,7 @@ subroutine fgenerate_fchl_acsf_and_gradients(coordinates, nuclear_charges, eleme & + angular * radial(l) * atm * rdecay(i,j) * rdecay(i,k) do t = 1, 3 - + ! Add up all gradient contributions wrt atom i atom_grad(z:z + nabasis - 1, i, t) = atom_grad(z:z + nabasis - 1, i, t) + & & d_angular * d_angular_d_i(t) * radial(l) * atm * rdecay(i,j) * rdecay(i,k) + & @@ -1293,40 +1289,40 @@ subroutine fgenerate_fchl_acsf_and_gradients(coordinates, nuclear_charges, eleme & angular * radial(l) * (atm_i * d_atm_ki(t) + atm_j * d_atm_kj(t) & & + atm_k * d_atm_kk(t) + d_atm_extra_k(t)) * three_body_weight * rdecay(i,j) * rdecay(i,k) - & & angular * radial(l) * rdecay(i,j) * d_ikdecay(t) * atm - + enddo enddo enddo enddo rep(i, twobody_size + 1:) = rep(i, twobody_size + 1:) + atom_rep grad(i, twobody_size + 1:,:,:) = grad(i, twobody_size + 1:,:,:) + atom_grad - enddo +enddo !$OMP END PARALLEL DO - deallocate(rdecay) - deallocate(element_types) - deallocate(distance_matrix) - deallocate(inv_distance_matrix) - deallocate(sq_distance_matrix) - deallocate(inv_sq_distance_matrix) - deallocate(atom_rep) - deallocate(atom_grad) - deallocate(a) - deallocate(b) - deallocate(c) - deallocate(radial) - deallocate(angular_base) - deallocate(angular) - deallocate(d_angular) - deallocate(d_angular_d_i) - deallocate(d_angular_d_j) - deallocate(d_angular_d_k) - deallocate(d_radial) - deallocate(d_radial_d_i) - deallocate(d_radial_d_j) - deallocate(d_radial_d_k) - deallocate(d_ijdecay) - deallocate(d_ikdecay) +deallocate(rdecay) +deallocate(element_types) +deallocate(distance_matrix) +deallocate(inv_distance_matrix) +deallocate(sq_distance_matrix) +deallocate(inv_sq_distance_matrix) +deallocate(atom_rep) +deallocate(atom_grad) +deallocate(a) +deallocate(b) +deallocate(c) +deallocate(radial) +deallocate(angular_base) +deallocate(angular) +deallocate(d_angular) +deallocate(d_angular_d_i) +deallocate(d_angular_d_j) +deallocate(d_angular_d_k) +deallocate(d_radial) +deallocate(d_radial_d_i) +deallocate(d_radial_d_j) +deallocate(d_radial_d_k) +deallocate(d_ijdecay) +deallocate(d_ikdecay) end subroutine fgenerate_fchl_acsf_and_gradients diff --git a/qml/representations/representations.py b/qml/representations/representations.py index f3f9be9ed..35610d68b 100644 --- a/qml/representations/representations.py +++ b/qml/representations/representations.py @@ -685,7 +685,7 @@ def generate_fchl_acsf(nuclear_charges, coordinates, elements = [1,6,7,8,16], natoms_tot=natoms if cell is not None: - nExtend = (np.floor(rcut/np.linalg.norm(cell,2,axis = 0)) + 1).astype(int) + nExtend = (np.floor(max(rcut,acut)/np.linalg.norm(cell,2,axis = 0)) + 1).astype(int) true_coords=coordinates for i in range(-nExtend[0],nExtend[0] + 1): for j in range(-nExtend[1],nExtend[1] + 1): diff --git a/test/test_fchl_acsf_pbc.py b/test/test_fchl_acsf_pbc.py index 52c5ff102..77bf75c0e 100644 --- a/test/test_fchl_acsf_pbc.py +++ b/test/test_fchl_acsf_pbc.py @@ -35,6 +35,7 @@ REP_PARAMS = dict() REP_PARAMS["elements"] = [1, 6, 7, 8, 16] REP_PARAMS["rcut"]=5.0 +REP_PARAMS["acut"]=5.0 random.seed(1) cell_added_cutoff=0.1 @@ -55,7 +56,7 @@ def generate_fchl_acsf_brute_pbc(nuclear_charges, coordinates, cell, gradients=F num_atoms=len(nuclear_charges) all_coords=deepcopy(coordinates) all_charges=deepcopy(nuclear_charges) - nExtend = (np.floor(REP_PARAMS["rcut"]/np.linalg.norm(cell,2,axis = 0)) + 1).astype(int) + nExtend = (np.floor(max(REP_PARAMS["rcut"], REP_PARAMS["acut"])/np.linalg.norm(cell,2,axis = 0)) + 1).astype(int) print("Checked nExtend:", nExtend, ", gradient calculation:", gradients) for i in range(-nExtend[0],nExtend[0] + 1): for j in range(-nExtend[1],nExtend[1] + 1): @@ -83,7 +84,7 @@ def test_fchl_acsf_pbc(): qm7_dir = os.path.dirname(os.path.realpath(__file__))+"/qm7" os.chdir(qm7_dir) all_xyzs=os.listdir() - test_xyzs=random.sample(all_xyzs, 3) + test_xyzs=random.sample(all_xyzs, 5) reps_no_grad1=[] reps_no_grad2=[] From 8bcc070959713da902bcb59f91aae3e872964ef0 Mon Sep 17 00:00:00 2001 From: kvkarandashev Date: Mon, 1 Apr 2024 16:28:38 +0200 Subject: [PATCH 6/8] Minor nitpick. --- qml/representations/facsf.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/qml/representations/facsf.f90 b/qml/representations/facsf.f90 index a918b1793..dc51b3d01 100644 --- a/qml/representations/facsf.f90 +++ b/qml/representations/facsf.f90 @@ -1063,8 +1063,8 @@ subroutine fgenerate_fchl_acsf_and_gradients(coordinates, nuclear_charges, eleme !$OMP PARALLEL DO PRIVATE(m) SCHEDULE(dynamic) do j = 2, natoms do i = 1, j-1 - m = element_types(i) if (distance_matrix(i,j)>rcut) cycle + m = element_types(i) rep(j, (m-1)*nbasis2 + 1:m*nbasis2) = rep(j, (m-1)*nbasis2 + 1:m*nbasis2) + add_rep(:, i, j) do k=1, 3 grad(j, (m-1)*nbasis2 + 1:m*nbasis2, j, k) = grad(j, (m-1)*nbasis2 + 1:m*nbasis2, j, k) - add_grad(:, k, i, j) From 588e96731354ef6679e037610030baac4baa6307 Mon Sep 17 00:00:00 2001 From: kvkarandashev Date: Thu, 11 Apr 2024 15:21:09 +0200 Subject: [PATCH 7/8] PBC gradient fix. --- qml/representations/facsf.f90 | 57 +++++++++++++++++++++-------------- test/test_fchl_acsf_pbc.py | 19 ++++++++++-- 2 files changed, 52 insertions(+), 24 deletions(-) diff --git a/qml/representations/facsf.f90 b/qml/representations/facsf.f90 index dc51b3d01..2a9784881 100644 --- a/qml/representations/facsf.f90 +++ b/qml/representations/facsf.f90 @@ -96,6 +96,17 @@ function calc_cos_angle(a, b, c) result(cos_angle) end function calc_cos_angle +function true_atom_id(atom_id, natoms) result(tatom_id) + integer, intent(in):: atom_id, natoms + integer:: tatom_id + + if (atom_id <= tatom_id) then + tatom_id=atom_id + else + tatom_id=modulo(atom_id-1, natoms)+1 + endif +end function true_atom_id + end module acsf_utils @@ -637,7 +648,7 @@ subroutine fgenerate_fchl_acsf(coordinates, nuclear_charges, elements, & & Rs2, Rs3, Ts, eta2, eta3, zeta, rcut, acut, natoms, natoms_tot, rep_size, & & two_body_decay, three_body_decay, three_body_weight, rep) - use acsf_utils, only: decay, calc_angle, calc_cos_angle + use acsf_utils, only: decay, calc_angle, calc_cos_angle, true_atom_id implicit none @@ -661,7 +672,7 @@ subroutine fgenerate_fchl_acsf(coordinates, nuclear_charges, elements, & double precision, intent(out), dimension(natoms, rep_size) :: rep integer:: two_body_rep_size - integer :: i, j, k, l, n, m, o, p, q, s, z, nelements, nbasis2, nbasis3, nabasis, j_id, j_id_max + integer :: i, j, k, l, n, m, o, p, q, s, z, nelements, nbasis2, nbasis3, nabasis, j_id, j_id_max, true_i integer, allocatable, dimension(:) :: element_types double precision :: rij, rik, angle, cos_1, cos_2, cos_3, invcut ! double precision :: angle_1, angle_2, angle_3 @@ -691,8 +702,9 @@ subroutine fgenerate_fchl_acsf(coordinates, nuclear_charges, elements, & ! Store element index of every atom !$OMP PARALLEL DO SCHEDULE(dynamic) do i = 1, natoms_tot + true_i=true_atom_id(i, natoms) do j = 1, nelements - if (nuclear_charges(modulo(i-1, natoms)+1) .eq. elements(j)) then + if (nuclear_charges(true_i) .eq. elements(j)) then element_types(i) = j exit endif @@ -869,7 +881,7 @@ subroutine fgenerate_fchl_acsf_and_gradients(coordinates, nuclear_charges, eleme & Rs2, Rs3, Ts, eta2, eta3, zeta, rcut, acut, natoms, natoms_tot, rep_size, & & two_body_decay, three_body_decay, three_body_weight, rep, grad) - use acsf_utils, only: decay, calc_angle, calc_cos_angle + use acsf_utils, only: decay, calc_angle, calc_cos_angle, true_atom_id implicit none @@ -903,7 +915,7 @@ subroutine fgenerate_fchl_acsf_and_gradients(coordinates, nuclear_charges, eleme double precision, dimension(:, :, :), allocatable:: add_rep double precision, dimension(:, :, :, :), allocatable :: add_grad - integer :: i, j, k, l, m, n, p, q, s, t, z, nelements, nbasis2, nbasis3, nabasis, twobody_size + integer :: i, j, k, l, m, n, p, q, s, t, z, nelements, nbasis2, nbasis3, nabasis, twobody_size, true_j, true_k integer, allocatable, dimension(:) :: element_types double precision :: rij, rik, angle, dot, rij2, rik2, invrij, invrik, invrij2, invrik2, invcut double precision, allocatable, dimension(:) :: radial_base, radial, angular, a, b, c, atom_rep @@ -937,13 +949,13 @@ subroutine fgenerate_fchl_acsf_and_gradients(coordinates, nuclear_charges, eleme ! Number of unique elements nelements = size(elements) ! Allocate temporary - allocate(element_types(natoms_tot)) + allocate(element_types(natoms)) ! Store element index of every atom !$OMP PARALLEL DO SCHEDULE(dynamic) - do i = 1, natoms_tot + do i = 1, natoms do j = 1, nelements - if (nuclear_charges(modulo(i-1, natoms)+1) .eq. elements(j)) then + if (nuclear_charges(i) .eq. elements(j)) then element_types(i) = j exit endif @@ -1013,12 +1025,14 @@ subroutine fgenerate_fchl_acsf_and_gradients(coordinates, nuclear_charges, eleme !$OMP scaling, radial_base, radial, dx, part, dscal, ddecay) SCHEDULE(dynamic) do i = 1, natoms ! The element index of atom i - m = element_types(i) ! moved inside loop to enable COLLAPSE + m = element_types(i) do j = i + 1, natoms_tot rij = distance_matrix(i,j) if (rij > rcut) continue + true_j=true_atom_id(i, natoms) + if (true_j < i) continue ! The element index of atom j - n = element_types(j) + n = element_types(true_j) ! Distance between atoms i and j invrij = inv_distance_matrix(i,j) mu = log(rij / sqrt(1.0d0 + eta2 * inv_sq_distance_matrix(i, j))) @@ -1034,7 +1048,7 @@ subroutine fgenerate_fchl_acsf_and_gradients(coordinates, nuclear_charges, eleme radial(:) = radial_base(:) * scaling * rdecay(i,j) rep(i, (n-1)*nbasis2 + 1:n*nbasis2) = rep(i, (n-1)*nbasis2 + 1:n*nbasis2) + radial - if (j<=natoms) add_rep(:, i, j)=radial + add_rep(:, i, true_j)=add_rep(:, i, true_j)+radial do k = 1, 3 dx = -(coordinates(i,k) - coordinates(j,k)) part(:) = ((log_Rs2(:) - mu) * (-dx *(rij**2 * exp_s2 + eta2) / (rij * sqrt(exp_s2))**3) & @@ -1050,11 +1064,9 @@ subroutine fgenerate_fchl_acsf_and_gradients(coordinates, nuclear_charges, eleme ! The gradients wrt coordinates grad(i, (n-1)*nbasis2 + 1:n*nbasis2, i, k) = grad(i, (n-1)*nbasis2 + 1:n*nbasis2, i, k) + part - if (j<=natoms) then - grad(i, (n-1)*nbasis2 + 1:n*nbasis2, j, k) = grad(i, (n-1)*nbasis2 + 1:n*nbasis2, j, k) - part - grad(j, (m-1)*nbasis2 + 1:m*nbasis2, i, k) = grad(j, (m-1)*nbasis2 + 1:m*nbasis2, i, k) + part - add_grad(:, k, i, j)=part - endif + grad(i, (n-1)*nbasis2 + 1:n*nbasis2, true_j, k) = grad(i, (n-1)*nbasis2 + 1:n*nbasis2, true_j, k) - part + grad(true_j, (m-1)*nbasis2 + 1:m*nbasis2, i, k) = grad(true_j, (m-1)*nbasis2 + 1:m*nbasis2, i, k) + part + add_grad(:, k, i, true_j)=add_grad(:, k, i, true_j)+part enddo enddo enddo @@ -1152,7 +1164,8 @@ subroutine fgenerate_fchl_acsf_and_gradients(coordinates, nuclear_charges, eleme rij = distance_matrix(i,j) if (rij > acut) cycle ! index of the element of atom j - n = element_types(j) + true_j=true_atom_id(j, natoms) + n = element_types(true_j) ! squared distance between atom i and j rij2 = sq_distance_matrix(i,j) ! inverse distance between atom i and j @@ -1165,7 +1178,8 @@ subroutine fgenerate_fchl_acsf_and_gradients(coordinates, nuclear_charges, eleme rik = distance_matrix(i,k) if (rik > acut) cycle ! index of the element of atom k - m = element_types(k) + true_k=true_atom_id(k, natoms) + m = element_types(true_k) ! squared distance between atom i and k rik2 = sq_distance_matrix(i,k) ! inverse distance between atom i and k @@ -1254,6 +1268,7 @@ subroutine fgenerate_fchl_acsf_and_gradients(coordinates, nuclear_charges, eleme ! Get index of where the contributions of atoms i,j,k should be added s = nbasis3 * nabasis * (-(p * (p + 1))/2 + q + nelements * p) + 1 + do l = 1, nbasis3 ! Get index of where the contributions of atoms i,j,k should be added @@ -1273,8 +1288,7 @@ subroutine fgenerate_fchl_acsf_and_gradients(coordinates, nuclear_charges, eleme & + atm_k * d_atm_ik(t) + d_atm_extra_i(t)) * three_body_weight * rdecay(i,j) * rdecay(i,k) + & & angular * radial(l) * (d_ijdecay(t) * rdecay(i,k) + rdecay(i,j) * d_ikdecay(t)) * atm ! Add up all gradient contributions wrt atom j - if (j<=natoms) & - atom_grad(z:z + nabasis - 1, j, t) = atom_grad(z:z + nabasis - 1, j, t) + & + atom_grad(z:z + nabasis - 1, true_j, t) = atom_grad(z:z + nabasis - 1, true_j, t) + & & d_angular * d_angular_d_j(t) * radial(l) * atm * rdecay(i,j) * rdecay(i,k) + & & angular * d_radial(l) * d_radial_d_j(t) * atm * rdecay(i,j) * rdecay(i,k) + & & angular * radial(l) * (atm_i * d_atm_ji(t) + atm_j * d_atm_jj(t) & @@ -1282,8 +1296,7 @@ subroutine fgenerate_fchl_acsf_and_gradients(coordinates, nuclear_charges, eleme & angular * radial(l) * d_ijdecay(t) * rdecay(i,k) * atm ! Add up all gradient contributions wrt atom k - if (k<=natoms) & - atom_grad(z:z + nabasis - 1, k, t) = atom_grad(z:z + nabasis - 1, k, t) + & + atom_grad(z:z + nabasis - 1, true_k, t) = atom_grad(z:z + nabasis - 1, true_k, t) + & & d_angular * d_angular_d_k(t) * radial(l) * atm * rdecay(i,j) * rdecay(i,k) + & & angular * d_radial(l) * d_radial_d_k(t) * atm * rdecay(i,j) * rdecay(i,k) + & & angular * radial(l) * (atm_i * d_atm_ki(t) + atm_j * d_atm_kj(t) & diff --git a/test/test_fchl_acsf_pbc.py b/test/test_fchl_acsf_pbc.py index 77bf75c0e..64c52b983 100644 --- a/test/test_fchl_acsf_pbc.py +++ b/test/test_fchl_acsf_pbc.py @@ -52,6 +52,19 @@ def suitable_cell(coords): min_coords=np.minimum(min_coords, atom_coords) return np.diag((max_coords-min_coords)*(1.0+cell_added_cutoff)) +def pbc_corrected_drep(drep, num_atoms): + new_shape=list(drep.shape) + new_shape[0]=num_atoms + new_shape[2]=num_atoms + new_drep=np.zeros(new_shape) + num_atoms_tot=drep.shape[0] + for i in range(num_atoms): + for j in range(num_atoms_tot): + true_j=j % num_atoms + new_drep[i, :, true_j, :] +=drep[i, :, j, :] + return new_drep + + def generate_fchl_acsf_brute_pbc(nuclear_charges, coordinates, cell, gradients=False): num_atoms=len(nuclear_charges) all_coords=deepcopy(coordinates) @@ -70,7 +83,7 @@ def generate_fchl_acsf_brute_pbc(nuclear_charges, coordinates, cell, gradients=F rep=generate_fchl_acsf(all_charges, all_coords, gradients=False, **REP_PARAMS) rep=rep[:num_atoms,:] if gradients: - drep=drep[:num_atoms,:,:num_atoms,:] + drep= pbc_corrected_drep(drep, num_atoms) return rep, drep else: return rep @@ -79,12 +92,14 @@ def ragged_array_close(arr1, arr2, error_msg): for el1, el2 in zip(arr1, arr2): assert np.allclose(el1, el2), error_msg + + def test_fchl_acsf_pbc(): qm7_dir = os.path.dirname(os.path.realpath(__file__))+"/qm7" os.chdir(qm7_dir) all_xyzs=os.listdir() - test_xyzs=random.sample(all_xyzs, 5) + test_xyzs=random.sample(all_xyzs, 16) reps_no_grad1=[] reps_no_grad2=[] From 7d13151629da7c580e3c509b414c6ed50d2cb382 Mon Sep 17 00:00:00 2001 From: kvkarandashev Date: Fri, 12 Apr 2024 15:37:54 +0200 Subject: [PATCH 8/8] PBC gradient fix. --- qml/representations/facsf.f90 | 2484 ++++++++++++++++----------------- 1 file changed, 1218 insertions(+), 1266 deletions(-) diff --git a/qml/representations/facsf.f90 b/qml/representations/facsf.f90 index 2a9784881..7d60efd0a 100644 --- a/qml/representations/facsf.f90 +++ b/qml/representations/facsf.f90 @@ -1,303 +1,278 @@ - - - - - - - - - - - - - - - - - - - - - module acsf_utils - implicit none + implicit none contains -function decay(r, invrc, natoms) result(f) + function decay(r, invrc, natoms) result(f) - implicit none + implicit none - double precision, intent(in), dimension(:,:) :: r - double precision, intent(in) :: invrc - integer, intent(in) :: natoms - double precision, dimension(natoms, natoms) :: f + double precision, intent(in), dimension(:, :) :: r + double precision, intent(in) :: invrc + integer, intent(in) :: natoms + double precision, dimension(natoms, natoms) :: f - double precision, parameter :: pi = 4.0d0 * atan(1.0d0) + double precision, parameter :: pi = 4.0d0*atan(1.0d0) - ! Decaying function reaching 0 at rc - f = 0.5d0 * (cos(pi * r * invrc) + 1.0d0) + ! Decaying function reaching 0 at rc + f = 0.5d0*(cos(pi*r*invrc) + 1.0d0) + end function decay -end function decay + function calc_angle(a, b, c) result(angle) -function calc_angle(a, b, c) result(angle) + implicit none - implicit none + double precision, intent(in), dimension(3) :: a + double precision, intent(in), dimension(3) :: b + double precision, intent(in), dimension(3) :: c - double precision, intent(in), dimension(3) :: a - double precision, intent(in), dimension(3) :: b - double precision, intent(in), dimension(3) :: c + double precision, dimension(3) :: v1 + double precision, dimension(3) :: v2 - double precision, dimension(3) :: v1 - double precision, dimension(3) :: v2 + double precision :: cos_angle + double precision :: angle - double precision :: cos_angle - double precision :: angle + v1 = a - b + v2 = c - b - v1 = a - b - v2 = c - b + v1 = v1/norm2(v1) + v2 = v2/norm2(v2) - v1 = v1 / norm2(v1) - v2 = v2 / norm2(v2) + cos_angle = dot_product(v1, v2) - cos_angle = dot_product(v1,v2) + ! Clipping + if (cos_angle > 1.0d0) cos_angle = 1.0d0 + if (cos_angle < -1.0d0) cos_angle = -1.0d0 - ! Clipping - if (cos_angle > 1.0d0) cos_angle = 1.0d0 - if (cos_angle < -1.0d0) cos_angle = -1.0d0 + angle = acos(cos_angle) - angle = acos(cos_angle) + end function calc_angle -end function calc_angle + function calc_cos_angle(a, b, c) result(cos_angle) -function calc_cos_angle(a, b, c) result(cos_angle) + implicit none - implicit none + double precision, intent(in), dimension(3) :: a + double precision, intent(in), dimension(3) :: b + double precision, intent(in), dimension(3) :: c - double precision, intent(in), dimension(3) :: a - double precision, intent(in), dimension(3) :: b - double precision, intent(in), dimension(3) :: c + double precision, dimension(3) :: v1 + double precision, dimension(3) :: v2 - double precision, dimension(3) :: v1 - double precision, dimension(3) :: v2 + double precision :: cos_angle - double precision :: cos_angle + v1 = a - b + v2 = c - b - v1 = a - b - v2 = c - b + v1 = v1/norm2(v1) + v2 = v2/norm2(v2) - v1 = v1 / norm2(v1) - v2 = v2 / norm2(v2) + cos_angle = dot_product(v1, v2) - cos_angle = dot_product(v1,v2) + end function calc_cos_angle -end function calc_cos_angle + function true_atom_id(atom_id, natoms) result(tatom_id) + integer, intent(in):: atom_id, natoms + integer:: tatom_id -function true_atom_id(atom_id, natoms) result(tatom_id) - integer, intent(in):: atom_id, natoms - integer:: tatom_id + if (atom_id <= natoms) then + tatom_id = atom_id + else + tatom_id = modulo(atom_id - 1, natoms) + 1 + end if - if (atom_id <= tatom_id) then - tatom_id=atom_id - else - tatom_id=modulo(atom_id-1, natoms)+1 - endif -end function true_atom_id + end function end module acsf_utils - subroutine fgenerate_acsf(coordinates, nuclear_charges, elements, & & Rs2, Rs3, Ts, eta2, eta3, zeta, rcut, acut, natoms, rep_size, rep) - use acsf_utils, only: decay, calc_angle - - implicit none - - double precision, intent(in), dimension(:, :) :: coordinates - integer, intent(in), dimension(:) :: nuclear_charges - integer, intent(in), dimension(:) :: elements - double precision, intent(in), dimension(:) :: Rs2 - double precision, intent(in), dimension(:) :: Rs3 - double precision, intent(in), dimension(:) :: Ts - double precision, intent(in) :: eta2 - double precision, intent(in) :: eta3 - double precision, intent(in) :: zeta - double precision, intent(in) :: rcut - double precision, intent(in) :: acut - integer, intent(in) :: natoms - integer, intent(in) :: rep_size - double precision, intent(out), dimension(natoms, rep_size) :: rep - - integer :: i, j, k, l, n, m, p, q, s, z, nelements, nbasis2, nbasis3, nabasis - integer, allocatable, dimension(:) :: element_types - double precision :: rij, rik, angle, invcut - double precision, allocatable, dimension(:) :: radial, angular, a, b, c - double precision, allocatable, dimension(:, :) :: distance_matrix, rdecay, rep_subset - - double precision, parameter :: pi = 4.0d0 * atan(1.0d0) - - if (natoms /= size(nuclear_charges, dim=1)) then - write(*,*) "ERROR: Atom Centered Symmetry Functions creation" - write(*,*) natoms, "coordinates, but", & - & size(nuclear_charges, dim=1), "atom_types!" - stop - endif - - - ! number of element types - nelements = size(elements) - ! Allocate temporary - allocate(element_types(natoms)) - - ! Store element index of every atom - !$OMP PARALLEL DO - do i = 1, natoms - do j = 1, nelements - if (nuclear_charges(i) .eq. elements(j)) then - element_types(i) = j - cycle - endif - enddo - enddo - !$OMP END PARALLEL DO - - - ! Get distance matrix - ! Allocate temporary - allocate(distance_matrix(natoms, natoms)) - distance_matrix = 0.0d0 - - - !$OMP PARALLEL DO PRIVATE(rij) - do i = 1, natoms - do j = i+1, natoms - rij = norm2(coordinates(j,:) - coordinates(i,:)) - distance_matrix(i, j) = rij - distance_matrix(j, i) = rij - enddo - enddo - !$OMP END PARALLEL DO - - ! number of basis functions in the two body term - nbasis2 = size(Rs2) - - ! Inverse of the two body cutoff - invcut = 1.0d0 / rcut - ! pre-calculate the radial decay in the two body terms - rdecay = decay(distance_matrix, invcut, natoms) - - ! Allocate temporary - allocate(radial(nbasis2)) - allocate(rep_subset(natoms, nelements * nbasis2)) - - rep_subset = 0.0d0 - - !$OMP PARALLEL DO PRIVATE(n,m,rij,radial) COLLAPSE(2) REDUCTION(+:rep_subset) SCHEDULE(dynamic) - do i = 1, natoms - do j = 1, natoms - if (j .le. i) cycle - rij = distance_matrix(i,j) - if (rij <= rcut) then - ! index of the element of atom i - m = element_types(i) - ! index of the element of atom j - n = element_types(j) - ! distance between atoms i and j - ! two body term of the representation - radial = exp(-eta2*(rij - Rs2)**2) * rdecay(i,j) - rep_subset(i, (n-1)*nbasis2 + 1:n*nbasis2) = rep_subset(i, (n-1)*nbasis2 + 1:n*nbasis2) + radial - rep_subset(j, (m-1)*nbasis2 + 1:m*nbasis2) = rep_subset(j, (m-1)*nbasis2 + 1:m*nbasis2) + radial - endif - enddo - enddo - !$OMP END PARALLEL DO - - rep(:, 1:nelements * nbasis2) = rep_subset(:,:) - - deallocate(radial) - deallocate(rep_subset) - - ! number of radial basis functions in the three body term - nbasis3 = size(Rs3) - ! number of radial basis functions in the three body term - nabasis = size(Ts) - - ! Inverse of the three body cutoff - invcut = 1.0d0 / acut - ! pre-calculate the radial decay in the three body terms - rdecay = decay(distance_matrix, invcut, natoms) - - ! Allocate temporary - allocate(rep_subset(natoms,rep_size - nelements * nbasis2)) - allocate(a(3)) - allocate(b(3)) - allocate(c(3)) - allocate(radial(nbasis3)) - allocate(angular(nabasis)) - - rep_subset = 0.0d0 - - ! This could probably be done more efficiently if it's a bottleneck - ! Also the order is a bit wobbly compared to the tensorflow implementation - !$OMP PARALLEL DO PRIVATE(rij, n, rik, m, a, b, c, angle, radial, angular, & - !$OMP p, q, s, z) REDUCTION(+:rep_subset) COLLAPSE(2) SCHEDULE(dynamic) - do i = 1, natoms - do j = 1, natoms - 1 - if (i .eq. j) cycle - ! distance between atoms i and j - rij = distance_matrix(i,j) - if (rij > acut) cycle + use acsf_utils, only: decay, calc_angle + + implicit none + + double precision, intent(in), dimension(:, :) :: coordinates + integer, intent(in), dimension(:) :: nuclear_charges + integer, intent(in), dimension(:) :: elements + double precision, intent(in), dimension(:) :: Rs2 + double precision, intent(in), dimension(:) :: Rs3 + double precision, intent(in), dimension(:) :: Ts + double precision, intent(in) :: eta2 + double precision, intent(in) :: eta3 + double precision, intent(in) :: zeta + double precision, intent(in) :: rcut + double precision, intent(in) :: acut + integer, intent(in) :: natoms + integer, intent(in) :: rep_size + double precision, intent(out), dimension(natoms, rep_size) :: rep + + integer :: i, j, k, l, n, m, p, q, s, z, nelements, nbasis2, nbasis3, nabasis + integer, allocatable, dimension(:) :: element_types + double precision :: rij, rik, angle, invcut + double precision, allocatable, dimension(:) :: radial, angular, a, b, c + double precision, allocatable, dimension(:, :) :: distance_matrix, rdecay, rep_subset + + double precision, parameter :: pi = 4.0d0*atan(1.0d0) + + if (natoms /= size(nuclear_charges, dim=1)) then + write (*, *) "ERROR: Atom Centered Symmetry Functions creation" + write (*, *) natoms, "coordinates, but", & + & size(nuclear_charges, dim=1), "atom_types!" + stop + end if + + ! number of element types + nelements = size(elements) + ! Allocate temporary + allocate (element_types(natoms)) + + ! Store element index of every atom + !$OMP PARALLEL DO + do i = 1, natoms + do j = 1, nelements + if (nuclear_charges(i) .eq. elements(j)) then + element_types(i) = j + cycle + end if + end do + end do + !$OMP END PARALLEL DO + + ! Get distance matrix + ! Allocate temporary + allocate (distance_matrix(natoms, natoms)) + distance_matrix = 0.0d0 + + !$OMP PARALLEL DO PRIVATE(rij) + do i = 1, natoms + do j = i + 1, natoms + rij = norm2(coordinates(j, :) - coordinates(i, :)) + distance_matrix(i, j) = rij + distance_matrix(j, i) = rij + end do + end do + !$OMP END PARALLEL DO + + ! number of basis functions in the two body term + nbasis2 = size(Rs2) + + ! Inverse of the two body cutoff + invcut = 1.0d0/rcut + ! pre-calculate the radial decay in the two body terms + rdecay = decay(distance_matrix, invcut, natoms) + + ! Allocate temporary + allocate (radial(nbasis2)) + allocate (rep_subset(natoms, nelements*nbasis2)) + + rep_subset = 0.0d0 + + !$OMP PARALLEL DO PRIVATE(n,m,rij,radial) COLLAPSE(2) REDUCTION(+:rep_subset) SCHEDULE(dynamic) + do i = 1, natoms + do j = 1, natoms + if (j .le. i) cycle + rij = distance_matrix(i, j) + if (rij <= rcut) then + ! index of the element of atom i + m = element_types(i) ! index of the element of atom j n = element_types(j) - do k = j + 1, natoms - if (i .eq. k) cycle - ! distance between atoms i and k - rik = distance_matrix(i,k) - if (rik > acut) cycle - ! index of the element of atom k - m = element_types(k) - ! coordinates of atoms j, i, k - a = coordinates(j,:) - b = coordinates(i,:) - c = coordinates(k,:) - ! angle between atoms i, j and k centered on i - angle = calc_angle(a,b,c) - ! The radial part of the three body terms including decay - radial = exp(-eta3*(0.5d0 * (rij+rik) - Rs3)**2) * rdecay(i,j) * rdecay(i,k) - ! The angular part of the three body terms - angular = 2.0d0 * ((1.0d0 + cos(angle - Ts)) * 0.5d0) ** zeta - ! The lowest of the element indices for atoms j and k - p = min(n,m) - 1 - ! The highest of the element indices for atoms j and k - q = max(n,m) - 1 - ! calculate the indices that the three body terms should be added to - s = nbasis3 * nabasis * (-(p * (p + 1))/2 + q + nelements * p) + 1 - do l = 1, nbasis3 - ! calculate the indices that the three body terms should be added to - z = s + (l-1) * nabasis - ! Add the contributions from atoms i,j and k - rep_subset(i, z:z + nabasis - 1) = rep_subset(i, z:z + nabasis - 1) + angular * radial(l) - enddo - enddo - enddo - enddo - !$OMP END PARALLEL DO - - rep(:, nelements * nbasis2 + 1:) = rep_subset(:,:) - - deallocate(element_types) - deallocate(rdecay) - deallocate(distance_matrix) - deallocate(rep_subset) - deallocate(a) - deallocate(b) - deallocate(c) - deallocate(radial) - deallocate(angular) + ! distance between atoms i and j + ! two body term of the representation + radial = exp(-eta2*(rij - Rs2)**2)*rdecay(i, j) + rep_subset(i, (n - 1)*nbasis2 + 1:n*nbasis2) = rep_subset(i, (n - 1)*nbasis2 + 1:n*nbasis2) + radial + rep_subset(j, (m - 1)*nbasis2 + 1:m*nbasis2) = rep_subset(j, (m - 1)*nbasis2 + 1:m*nbasis2) + radial + end if + end do + end do + !$OMP END PARALLEL DO + + rep(:, 1:nelements*nbasis2) = rep_subset(:, :) + + deallocate (radial) + deallocate (rep_subset) + + ! number of radial basis functions in the three body term + nbasis3 = size(Rs3) + ! number of radial basis functions in the three body term + nabasis = size(Ts) + + ! Inverse of the three body cutoff + invcut = 1.0d0/acut + ! pre-calculate the radial decay in the three body terms + rdecay = decay(distance_matrix, invcut, natoms) + + ! Allocate temporary + allocate (rep_subset(natoms, rep_size - nelements*nbasis2)) + allocate (a(3)) + allocate (b(3)) + allocate (c(3)) + allocate (radial(nbasis3)) + allocate (angular(nabasis)) + + rep_subset = 0.0d0 + + ! This could probably be done more efficiently if it's a bottleneck + ! Also the order is a bit wobbly compared to the tensorflow implementation + !$OMP PARALLEL DO PRIVATE(rij, n, rik, m, a, b, c, angle, radial, angular, & + !$OMP p, q, s, z) REDUCTION(+:rep_subset) COLLAPSE(2) SCHEDULE(dynamic) + do i = 1, natoms + do j = 1, natoms - 1 + if (i .eq. j) cycle + ! distance between atoms i and j + rij = distance_matrix(i, j) + if (rij > acut) cycle + ! index of the element of atom j + n = element_types(j) + do k = j + 1, natoms + if (i .eq. k) cycle + ! distance between atoms i and k + rik = distance_matrix(i, k) + if (rik > acut) cycle + ! index of the element of atom k + m = element_types(k) + ! coordinates of atoms j, i, k + a = coordinates(j, :) + b = coordinates(i, :) + c = coordinates(k, :) + ! angle between atoms i, j and k centered on i + angle = calc_angle(a, b, c) + ! The radial part of the three body terms including decay + radial = exp(-eta3*(0.5d0*(rij + rik) - Rs3)**2)*rdecay(i, j)*rdecay(i, k) + ! The angular part of the three body terms + angular = 2.0d0*((1.0d0 + cos(angle - Ts))*0.5d0)**zeta + ! The lowest of the element indices for atoms j and k + p = min(n, m) - 1 + ! The highest of the element indices for atoms j and k + q = max(n, m) - 1 + ! calculate the indices that the three body terms should be added to + s = nbasis3*nabasis*(-(p*(p + 1))/2 + q + nelements*p) + 1 + do l = 1, nbasis3 + ! calculate the indices that the three body terms should be added to + z = s + (l - 1)*nabasis + ! Add the contributions from atoms i,j and k + rep_subset(i, z:z + nabasis - 1) = rep_subset(i, z:z + nabasis - 1) + angular*radial(l) + end do + end do + end do + end do + !$OMP END PARALLEL DO + + rep(:, nelements*nbasis2 + 1:) = rep_subset(:, :) + + deallocate (element_types) + deallocate (rdecay) + deallocate (distance_matrix) + deallocate (rep_subset) + deallocate (a) + deallocate (b) + deallocate (c) + deallocate (radial) + deallocate (angular) end subroutine fgenerate_acsf @@ -305,1037 +280,1014 @@ subroutine fgenerate_acsf_and_gradients(coordinates, nuclear_charges, elements, & Rs2, Rs3, Ts, eta2, eta3, zeta, rcut, acut, natoms, & & rep_size, rep, grad) - use acsf_utils, only: decay, calc_angle - - implicit none - - double precision, intent(in), dimension(:, :) :: coordinates - integer, intent(in), dimension(:) :: nuclear_charges - integer, intent(in), dimension(:) :: elements - double precision, intent(in), dimension(:) :: Rs2 - double precision, intent(in), dimension(:) :: Rs3 - double precision, intent(in), dimension(:) :: Ts - double precision, intent(in) :: eta2 - double precision, intent(in) :: eta3 - double precision, intent(in) :: zeta - double precision, intent(in) :: rcut - double precision, intent(in) :: acut - integer, intent(in) :: natoms - integer, intent(in) :: rep_size - double precision, intent(out), dimension(natoms, rep_size) :: rep - double precision, intent(out), dimension(natoms, rep_size, natoms, 3) :: grad - - integer :: i, j, k, l, n, m, p, q, s, t, z, nelements, nbasis2, nbasis3, nabasis, twobody_size - integer, allocatable, dimension(:) :: element_types - double precision :: rij, rik, angle, dot, rij2, rik2, invrij, invrik, invrij2, invrik2, invcut - double precision, allocatable, dimension(:) :: radial_base, radial, angular, a, b, c, atom_rep - double precision, allocatable, dimension(:) :: angular_base, d_radial, d_radial_d_i - double precision, allocatable, dimension(:) :: d_radial_d_j, d_radial_d_k, d_ijdecay, d_ikdecay - double precision, allocatable, dimension(:) :: d_angular, part, radial_part - double precision, allocatable, dimension(:) :: d_angular_d_i, d_angular_d_j, d_angular_d_k - double precision, allocatable, dimension(:, :) :: distance_matrix, rdecay, sq_distance_matrix - double precision, allocatable, dimension(:, :) :: inv_distance_matrix, inv_sq_distance_matrix - double precision, allocatable, dimension(:, :) :: rep_subset - double precision, allocatable, dimension(:, :, :) :: atom_grad - double precision, allocatable, dimension(:, :, :, :) :: grad_subset - - double precision, parameter :: pi = 4.0d0 * atan(1.0d0) - - if (natoms /= size(nuclear_charges, dim=1)) then - write(*,*) "ERROR: Atom Centered Symmetry Functions creation" - write(*,*) natoms, "coordinates, but", & - & size(nuclear_charges, dim=1), "atom_types!" - stop - endif - - - ! Number of unique elements - nelements = size(elements) - ! Allocate temporary - allocate(element_types(natoms)) - - ! Store element index of every atom - !$OMP PARALLEL DO - do i = 1, natoms - do j = 1, nelements - if (nuclear_charges(i) .eq. elements(j)) then - element_types(i) = j - cycle - endif - enddo - enddo - !$OMP END PARALLEL DO - - - - ! Get distance matrix - ! Allocate temporary - allocate(distance_matrix(natoms, natoms)) - allocate(sq_distance_matrix(natoms, natoms)) - allocate(inv_distance_matrix(natoms, natoms)) - allocate(inv_sq_distance_matrix(natoms, natoms)) - distance_matrix = 0.0d0 - sq_distance_matrix = 0.0d0 - inv_distance_matrix = 0.0d0 - inv_sq_distance_matrix = 0.0d0 - - - !$OMP PARALLEL DO PRIVATE(rij,rij2,invrij,invrij2) SCHEDULE(dynamic) - do i = 1, natoms - do j = 1, natoms - if (j .le. i) cycle - rij = norm2(coordinates(j,:) - coordinates(i,:)) - distance_matrix(i, j) = rij - distance_matrix(j, i) = rij - rij2 = rij * rij - sq_distance_matrix(i, j) = rij2 - sq_distance_matrix(j, i) = rij2 - invrij = 1.0d0 / rij - inv_distance_matrix(i, j) = invrij - inv_distance_matrix(j, i) = invrij - invrij2 = invrij * invrij - inv_sq_distance_matrix(i, j) = invrij2 - inv_sq_distance_matrix(j, i) = invrij2 - enddo - enddo - !$OMP END PARALLEL DO - - - ! Number of two body basis functions - nbasis2 = size(Rs2) - - ! Inverse of the two body cutoff distance - invcut = 1.0d0 / rcut - ! Pre-calculate the two body decay - rdecay = decay(distance_matrix, invcut, natoms) - - ! Allocate temporary - allocate(radial_base(nbasis2)) - allocate(radial(nbasis2)) - allocate(radial_part(nbasis2)) - allocate(part(nbasis2)) - allocate(rep_subset(natoms, nelements * nbasis2)) - allocate(grad_subset(natoms, nelements * nbasis2, natoms, 3)) - - grad_subset = 0.0d0 - rep_subset = 0.0d0 - - !$OMP PARALLEL DO PRIVATE(m,n,rij,invrij,radial_base,radial,radial_part,part) REDUCTION(+:rep_subset,grad_subset) & - !$OMP SCHEDULE(dynamic) - do i = 1, natoms - ! The element index of atom i - m = element_types(i) - do j = i + 1, natoms - ! The element index of atom j - n = element_types(j) - ! Distance between atoms i and j - rij = distance_matrix(i,j) - if (rij <= rcut) then - invrij = inv_distance_matrix(i,j) - ! part of the two body terms - radial_base = exp(-eta2*(rij - Rs2)**2) - ! The full two body term between atom i and j - radial = radial_base * rdecay(i,j) - ! Add the contributions from atoms i and j - rep_subset(i, (n-1)*nbasis2 + 1:n*nbasis2) = rep_subset(i, (n-1)*nbasis2 + 1:n*nbasis2) + radial - rep_subset(j, (m-1)*nbasis2 + 1:m*nbasis2) = rep_subset(j, (m-1)*nbasis2 + 1:m*nbasis2) + radial - ! Part of the gradients that can be reused for x,y and z coordinates - radial_part = - radial_base * invrij * (2.0d0 * eta2 * (rij - Rs2) * rdecay(i,j) + & - & 0.5d0 * pi * sin(pi*rij * invcut) * invcut) - do k = 1, 3 - ! The gradients wrt coordinates - part = radial_part * (coordinates(i,k) - coordinates(j,k)) - grad_subset(i, (n-1)*nbasis2 + 1:n*nbasis2, i, k) = & - grad_subset(i, (n-1)*nbasis2 + 1:n*nbasis2, i, k) + part - grad_subset(i, (n-1)*nbasis2 + 1:n*nbasis2, j, k) = & - grad_subset(i, (n-1)*nbasis2 + 1:n*nbasis2, j, k) - part - grad_subset(j, (m-1)*nbasis2 + 1:m*nbasis2, j, k) = & - grad_subset(j, (m-1)*nbasis2 + 1:m*nbasis2, j, k) - part - grad_subset(j, (m-1)*nbasis2 + 1:m*nbasis2, i, k) = & - grad_subset(j, (m-1)*nbasis2 + 1:m*nbasis2, i, k) + part - enddo - endif - enddo - enddo - !$OMP END PARALLEL DO - - rep(:,:nelements*nbasis2) = rep_subset(:,:) - grad(:,:nelements*nbasis2,:,:) = grad_subset(:,:,:,:) - - deallocate(radial_base) - deallocate(radial) - deallocate(radial_part) - deallocate(part) - deallocate(rep_subset) - deallocate(grad_subset) - - - ! Number of radial basis functions in the three body term - nbasis3 = size(Rs3) - ! Number of angular basis functions in the three body term - nabasis = size(Ts) - ! Size of two body terms - twobody_size = nelements * nbasis2 - - ! Inverse of the three body cutoff distance - invcut = 1.0d0 / acut - ! Pre-calculate the three body decay - rdecay = decay(distance_matrix, invcut, natoms) - - ! Allocate temporary - allocate(atom_rep(rep_size - twobody_size)) - allocate(atom_grad(rep_size - twobody_size, natoms, 3)) - allocate(a(3)) - allocate(b(3)) - allocate(c(3)) - allocate(radial(nbasis3)) - allocate(angular_base(nabasis)) - allocate(angular(nabasis)) - allocate(d_angular(nabasis)) - allocate(d_angular_d_i(3)) - allocate(d_angular_d_j(3)) - allocate(d_angular_d_k(3)) - allocate(d_radial(nbasis3)) - allocate(d_radial_d_i(3)) - allocate(d_radial_d_j(3)) - allocate(d_radial_d_k(3)) - allocate(d_ijdecay(3)) - allocate(d_ikdecay(3)) - - ! This could probably be done more efficiently if it's a bottleneck - ! The order is a bit wobbly compared to the tensorflow implementation - !$OMP PARALLEL DO PRIVATE(atom_rep,atom_grad,a,b,c,radial,angular_base, & - !$OMP angular,d_angular,rij,n,rij2,invrij,invrij2,d_angular_d_i, & - !$OMP d_angular_d_j,d_angular_d_k,rik,m,rik2,invrik,invrik2,angle, & - !$OMP p,q,dot,d_radial,d_radial_d_i,d_radial_d_j,d_radial_d_k,s,z, & - !$OMP d_ijdecay,d_ikdecay) SCHEDULE(dynamic) - do i = 1, natoms - atom_rep = 0.0d0 - atom_grad = 0.0d0 - do j = 1, natoms - 1 - if (i .eq. j) cycle - ! distance between atom i and j - rij = distance_matrix(i,j) - if (rij > acut) cycle - ! index of the element of atom j - n = element_types(j) - ! squared distance between atom i and j - rij2 = sq_distance_matrix(i,j) - ! inverse distance between atom i and j - invrij = inv_distance_matrix(i,j) - ! inverse squared distance between atom i and j - invrij2 = inv_sq_distance_matrix(i,j) - do k = j + 1, natoms - if (i .eq. k) cycle - ! distance between atom i and k - rik = distance_matrix(i,k) - if (rik > acut) cycle - ! index of the element of atom k - m = element_types(k) - ! squared distance between atom i and k - rik2 = sq_distance_matrix(i,k) - ! inverse distance between atom i and k - invrik = inv_distance_matrix(i,k) - ! inverse squared distance between atom i and k - invrik2 = inv_sq_distance_matrix(i,k) - ! coordinates of atoms j, i, k - a = coordinates(j,:) - b = coordinates(i,:) - c = coordinates(k,:) - ! angle between atom i, j and k, centered on i - angle = calc_angle(a,b,c) - ! part of the radial part of the 3body terms - radial = exp(-eta3*(0.5d0 * (rij+rik) - Rs3)**2) - ! used in the angular part of the 3body terms and in gradients - angular_base = ((1.0d0 + cos(angle - Ts)) * 0.5d0) - ! angular part of the 3body terms - angular = 2.0d0 * angular_base ** zeta - ! the lowest index of the elements of j,k - p = min(n,m) - 1 - ! the highest index of the elements of j,k - q = max(n,m) - 1 - ! Dot product between the vectors connecting atom i,j and i,k - dot = dot_product(a-b,c-b) - ! Part of the derivative of the angular basis functions wrt coordinates (dim(nabasis)) - ! including decay - d_angular = (zeta * angular * sin(angle-Ts) * rdecay(i,j) * rdecay(i,k)) / & - ! & (2.0d0 * sqrt(rij2 * rik2 - dot**2) * angular_base) - & (2.0d0 * max(1d-10, sqrt(abs(rij2 * rik2 - dot**2)) * angular_base)) - ! write(*,*) angular_base - ! Part of the derivative of the angular basis functions wrt atom j (dim(3)) - d_angular_d_j = c - b + dot * ((b - a) * invrij2) - ! Part of the derivative of the angular basis functions wrt atom k (dim(3)) - d_angular_d_k = a - b + dot * ((b - c) * invrik2) - ! Part of the derivative of the angular basis functions wrt atom i (dim(3)) - d_angular_d_i = - (d_angular_d_j + d_angular_d_k) - ! Part of the derivative of the radial basis functions wrt coordinates (dim(nbasis3)) - ! including decay - d_radial = radial * eta3 * (0.5d0 * (rij+rik) - Rs3) * rdecay(i,j) * rdecay(i,k) - ! Part of the derivative of the radial basis functions wrt atom j (dim(3)) - d_radial_d_j = (b - a) * invrij - ! Part of the derivative of the radial basis functions wrt atom k (dim(3)) - d_radial_d_k = (b - c) * invrik - ! Part of the derivative of the radial basis functions wrt atom i (dim(3)) - d_radial_d_i = - (d_radial_d_j + d_radial_d_k) - ! Part of the derivative of the i,j decay functions wrt coordinates (dim(3)) - d_ijdecay = - pi * (b - a) * sin(pi * rij * invcut) * 0.5d0 * invrij * invcut - ! Part of the derivative of the i,k decay functions wrt coordinates (dim(3)) - d_ikdecay = - pi * (b - c) * sin(pi * rik * invcut) * 0.5d0 * invrik * invcut - - ! Get index of where the contributions of atoms i,j,k should be added - s = nbasis3 * nabasis * (-(p * (p + 1))/2 + q + nelements * p) + 1 - do l = 1, nbasis3 - ! Get index of where the contributions of atoms i,j,k should be added - z = s + (l-1) * nabasis - ! Add the contributions for atoms i,j,k - atom_rep(z:z + nabasis - 1) = atom_rep(z:z + nabasis - 1) + angular * radial(l) * rdecay(i,j) * rdecay(i,k) - do t = 1, 3 - ! Add up all gradient contributions wrt atom i - atom_grad(z:z + nabasis - 1, i, t) = atom_grad(z:z + nabasis - 1, i, t) + & - & d_angular * d_angular_d_i(t) * radial(l) + & - & angular * d_radial(l) * d_radial_d_i(t) + & - & angular * radial(l) * (d_ijdecay(t) * rdecay(i,k) + rdecay(i,j) * d_ikdecay(t)) - ! Add up all gradient contributions wrt atom j - atom_grad(z:z + nabasis - 1, j, t) = atom_grad(z:z + nabasis - 1, j, t) + & - & d_angular * d_angular_d_j(t) * radial(l) + & - & angular * d_radial(l) * d_radial_d_j(t) - & - & angular * radial(l) * d_ijdecay(t) * rdecay(i,k) - ! Add up all gradient contributions wrt atom k - atom_grad(z:z + nabasis - 1, k, t) = atom_grad(z:z + nabasis - 1, k, t) + & - & d_angular * d_angular_d_k(t) * radial(l) + & - & angular * d_radial(l) * d_radial_d_k(t) - & - & angular * radial(l) * rdecay(i,j) * d_ikdecay(t) - enddo - enddo - enddo - enddo - rep(i, twobody_size + 1:) = atom_rep - grad(i, twobody_size + 1:,:,:) = atom_grad - enddo - !$OMP END PARALLEL DO - - - deallocate(rdecay) - deallocate(element_types) - deallocate(distance_matrix) - deallocate(inv_distance_matrix) - deallocate(sq_distance_matrix) - deallocate(inv_sq_distance_matrix) - deallocate(atom_rep) - deallocate(atom_grad) - deallocate(a) - deallocate(b) - deallocate(c) - deallocate(radial) - deallocate(angular_base) - deallocate(angular) - deallocate(d_angular) - deallocate(d_angular_d_i) - deallocate(d_angular_d_j) - deallocate(d_angular_d_k) - deallocate(d_radial) - deallocate(d_radial_d_i) - deallocate(d_radial_d_j) - deallocate(d_radial_d_k) - deallocate(d_ijdecay) - deallocate(d_ikdecay) - + use acsf_utils, only: decay, calc_angle + + implicit none + + double precision, intent(in), dimension(:, :) :: coordinates + integer, intent(in), dimension(:) :: nuclear_charges + integer, intent(in), dimension(:) :: elements + double precision, intent(in), dimension(:) :: Rs2 + double precision, intent(in), dimension(:) :: Rs3 + double precision, intent(in), dimension(:) :: Ts + double precision, intent(in) :: eta2 + double precision, intent(in) :: eta3 + double precision, intent(in) :: zeta + double precision, intent(in) :: rcut + double precision, intent(in) :: acut + integer, intent(in) :: natoms + integer, intent(in) :: rep_size + double precision, intent(out), dimension(natoms, rep_size) :: rep + double precision, intent(out), dimension(natoms, rep_size, natoms, 3) :: grad + + integer :: i, j, k, l, n, m, p, q, s, t, z, nelements, nbasis2, nbasis3, nabasis, twobody_size + integer, allocatable, dimension(:) :: element_types + double precision :: rij, rik, angle, dot, rij2, rik2, invrij, invrik, invrij2, invrik2, invcut + double precision, allocatable, dimension(:) :: radial_base, radial, angular, a, b, c, atom_rep + double precision, allocatable, dimension(:) :: angular_base, d_radial, d_radial_d_i + double precision, allocatable, dimension(:) :: d_radial_d_j, d_radial_d_k, d_ijdecay, d_ikdecay + double precision, allocatable, dimension(:) :: d_angular, part, radial_part + double precision, allocatable, dimension(:) :: d_angular_d_i, d_angular_d_j, d_angular_d_k + double precision, allocatable, dimension(:, :) :: distance_matrix, rdecay, sq_distance_matrix + double precision, allocatable, dimension(:, :) :: inv_distance_matrix, inv_sq_distance_matrix + double precision, allocatable, dimension(:, :) :: rep_subset + double precision, allocatable, dimension(:, :, :) :: atom_grad + double precision, allocatable, dimension(:, :, :, :) :: grad_subset + + double precision, parameter :: pi = 4.0d0*atan(1.0d0) + + if (natoms /= size(nuclear_charges, dim=1)) then + write (*, *) "ERROR: Atom Centered Symmetry Functions creation" + write (*, *) natoms, "coordinates, but", & + & size(nuclear_charges, dim=1), "atom_types!" + stop + end if + + ! Number of unique elements + nelements = size(elements) + ! Allocate temporary + allocate (element_types(natoms)) + + ! Store element index of every atom + !$OMP PARALLEL DO + do i = 1, natoms + do j = 1, nelements + if (nuclear_charges(i) .eq. elements(j)) then + element_types(i) = j + cycle + end if + end do + end do + !$OMP END PARALLEL DO + + ! Get distance matrix + ! Allocate temporary + allocate (distance_matrix(natoms, natoms)) + allocate (sq_distance_matrix(natoms, natoms)) + allocate (inv_distance_matrix(natoms, natoms)) + allocate (inv_sq_distance_matrix(natoms, natoms)) + distance_matrix = 0.0d0 + sq_distance_matrix = 0.0d0 + inv_distance_matrix = 0.0d0 + inv_sq_distance_matrix = 0.0d0 + + !$OMP PARALLEL DO PRIVATE(rij,rij2,invrij,invrij2) SCHEDULE(dynamic) + do i = 1, natoms + do j = 1, natoms + if (j .le. i) cycle + rij = norm2(coordinates(j, :) - coordinates(i, :)) + distance_matrix(i, j) = rij + distance_matrix(j, i) = rij + rij2 = rij*rij + sq_distance_matrix(i, j) = rij2 + sq_distance_matrix(j, i) = rij2 + invrij = 1.0d0/rij + inv_distance_matrix(i, j) = invrij + inv_distance_matrix(j, i) = invrij + invrij2 = invrij*invrij + inv_sq_distance_matrix(i, j) = invrij2 + inv_sq_distance_matrix(j, i) = invrij2 + end do + end do + !$OMP END PARALLEL DO + + ! Number of two body basis functions + nbasis2 = size(Rs2) + + ! Inverse of the two body cutoff distance + invcut = 1.0d0/rcut + ! Pre-calculate the two body decay + rdecay = decay(distance_matrix, invcut, natoms) + + ! Allocate temporary + allocate (radial_base(nbasis2)) + allocate (radial(nbasis2)) + allocate (radial_part(nbasis2)) + allocate (part(nbasis2)) + allocate (rep_subset(natoms, nelements*nbasis2)) + allocate (grad_subset(natoms, nelements*nbasis2, natoms, 3)) + + grad_subset = 0.0d0 + rep_subset = 0.0d0 + + !$OMP PARALLEL DO PRIVATE(m,n,rij,invrij,radial_base,radial,radial_part,part) REDUCTION(+:rep_subset,grad_subset) & + !$OMP SCHEDULE(dynamic) + do i = 1, natoms + ! The element index of atom i + m = element_types(i) + do j = i + 1, natoms + ! The element index of atom j + n = element_types(j) + ! Distance between atoms i and j + rij = distance_matrix(i, j) + if (rij <= rcut) then + invrij = inv_distance_matrix(i, j) + ! part of the two body terms + radial_base = exp(-eta2*(rij - Rs2)**2) + ! The full two body term between atom i and j + radial = radial_base*rdecay(i, j) + ! Add the contributions from atoms i and j + rep_subset(i, (n - 1)*nbasis2 + 1:n*nbasis2) = rep_subset(i, (n - 1)*nbasis2 + 1:n*nbasis2) + radial + rep_subset(j, (m - 1)*nbasis2 + 1:m*nbasis2) = rep_subset(j, (m - 1)*nbasis2 + 1:m*nbasis2) + radial + ! Part of the gradients that can be reused for x,y and z coordinates + radial_part = -radial_base*invrij*(2.0d0*eta2*(rij - Rs2)*rdecay(i, j) + & + & 0.5d0*pi*sin(pi*rij*invcut)*invcut) + do k = 1, 3 + ! The gradients wrt coordinates + part = radial_part*(coordinates(i, k) - coordinates(j, k)) + grad_subset(i, (n - 1)*nbasis2 + 1:n*nbasis2, i, k) = & + grad_subset(i, (n - 1)*nbasis2 + 1:n*nbasis2, i, k) + part + grad_subset(i, (n - 1)*nbasis2 + 1:n*nbasis2, j, k) = & + grad_subset(i, (n - 1)*nbasis2 + 1:n*nbasis2, j, k) - part + grad_subset(j, (m - 1)*nbasis2 + 1:m*nbasis2, j, k) = & + grad_subset(j, (m - 1)*nbasis2 + 1:m*nbasis2, j, k) - part + grad_subset(j, (m - 1)*nbasis2 + 1:m*nbasis2, i, k) = & + grad_subset(j, (m - 1)*nbasis2 + 1:m*nbasis2, i, k) + part + end do + end if + end do + end do + !$OMP END PARALLEL DO + + rep(:, :nelements*nbasis2) = rep_subset(:, :) + grad(:, :nelements*nbasis2, :, :) = grad_subset(:, :, :, :) + + deallocate (radial_base) + deallocate (radial) + deallocate (radial_part) + deallocate (part) + deallocate (rep_subset) + deallocate (grad_subset) + + ! Number of radial basis functions in the three body term + nbasis3 = size(Rs3) + ! Number of angular basis functions in the three body term + nabasis = size(Ts) + ! Size of two body terms + twobody_size = nelements*nbasis2 + + ! Inverse of the three body cutoff distance + invcut = 1.0d0/acut + ! Pre-calculate the three body decay + rdecay = decay(distance_matrix, invcut, natoms) + + ! Allocate temporary + allocate (atom_rep(rep_size - twobody_size)) + allocate (atom_grad(rep_size - twobody_size, natoms, 3)) + allocate (a(3)) + allocate (b(3)) + allocate (c(3)) + allocate (radial(nbasis3)) + allocate (angular_base(nabasis)) + allocate (angular(nabasis)) + allocate (d_angular(nabasis)) + allocate (d_angular_d_i(3)) + allocate (d_angular_d_j(3)) + allocate (d_angular_d_k(3)) + allocate (d_radial(nbasis3)) + allocate (d_radial_d_i(3)) + allocate (d_radial_d_j(3)) + allocate (d_radial_d_k(3)) + allocate (d_ijdecay(3)) + allocate (d_ikdecay(3)) + + ! This could probably be done more efficiently if it's a bottleneck + ! The order is a bit wobbly compared to the tensorflow implementation + !$OMP PARALLEL DO PRIVATE(atom_rep,atom_grad,a,b,c,radial,angular_base, & + !$OMP angular,d_angular,rij,n,rij2,invrij,invrij2,d_angular_d_i, & + !$OMP d_angular_d_j,d_angular_d_k,rik,m,rik2,invrik,invrik2,angle, & + !$OMP p,q,dot,d_radial,d_radial_d_i,d_radial_d_j,d_radial_d_k,s,z, & + !$OMP d_ijdecay,d_ikdecay) SCHEDULE(dynamic) + do i = 1, natoms + atom_rep = 0.0d0 + atom_grad = 0.0d0 + do j = 1, natoms - 1 + if (i .eq. j) cycle + ! distance between atom i and j + rij = distance_matrix(i, j) + if (rij > acut) cycle + ! index of the element of atom j + n = element_types(j) + ! squared distance between atom i and j + rij2 = sq_distance_matrix(i, j) + ! inverse distance between atom i and j + invrij = inv_distance_matrix(i, j) + ! inverse squared distance between atom i and j + invrij2 = inv_sq_distance_matrix(i, j) + do k = j + 1, natoms + if (i .eq. k) cycle + ! distance between atom i and k + rik = distance_matrix(i, k) + if (rik > acut) cycle + ! index of the element of atom k + m = element_types(k) + ! squared distance between atom i and k + rik2 = sq_distance_matrix(i, k) + ! inverse distance between atom i and k + invrik = inv_distance_matrix(i, k) + ! inverse squared distance between atom i and k + invrik2 = inv_sq_distance_matrix(i, k) + ! coordinates of atoms j, i, k + a = coordinates(j, :) + b = coordinates(i, :) + c = coordinates(k, :) + ! angle between atom i, j and k, centered on i + angle = calc_angle(a, b, c) + ! part of the radial part of the 3body terms + radial = exp(-eta3*(0.5d0*(rij + rik) - Rs3)**2) + ! used in the angular part of the 3body terms and in gradients + angular_base = ((1.0d0 + cos(angle - Ts))*0.5d0) + ! angular part of the 3body terms + angular = 2.0d0*angular_base**zeta + ! the lowest index of the elements of j,k + p = min(n, m) - 1 + ! the highest index of the elements of j,k + q = max(n, m) - 1 + ! Dot product between the vectors connecting atom i,j and i,k + dot = dot_product(a - b, c - b) + ! Part of the derivative of the angular basis functions wrt coordinates (dim(nabasis)) + ! including decay + d_angular = (zeta*angular*sin(angle - Ts)*rdecay(i, j)*rdecay(i, k))/ & + ! & (2.0d0 * sqrt(rij2 * rik2 - dot**2) * angular_base) + (2.0d0*max(1d-10, sqrt(abs(rij2*rik2 - dot**2))*angular_base)) + ! write(*,*) angular_base + ! Part of the derivative of the angular basis functions wrt atom j (dim(3)) + d_angular_d_j = c - b + dot*((b - a)*invrij2) + ! Part of the derivative of the angular basis functions wrt atom k (dim(3)) + d_angular_d_k = a - b + dot*((b - c)*invrik2) + ! Part of the derivative of the angular basis functions wrt atom i (dim(3)) + d_angular_d_i = -(d_angular_d_j + d_angular_d_k) + ! Part of the derivative of the radial basis functions wrt coordinates (dim(nbasis3)) + ! including decay + d_radial = radial*eta3*(0.5d0*(rij + rik) - Rs3)*rdecay(i, j)*rdecay(i, k) + ! Part of the derivative of the radial basis functions wrt atom j (dim(3)) + d_radial_d_j = (b - a)*invrij + ! Part of the derivative of the radial basis functions wrt atom k (dim(3)) + d_radial_d_k = (b - c)*invrik + ! Part of the derivative of the radial basis functions wrt atom i (dim(3)) + d_radial_d_i = -(d_radial_d_j + d_radial_d_k) + ! Part of the derivative of the i,j decay functions wrt coordinates (dim(3)) + d_ijdecay = -pi*(b - a)*sin(pi*rij*invcut)*0.5d0*invrij*invcut + ! Part of the derivative of the i,k decay functions wrt coordinates (dim(3)) + d_ikdecay = -pi*(b - c)*sin(pi*rik*invcut)*0.5d0*invrik*invcut + + ! Get index of where the contributions of atoms i,j,k should be added + s = nbasis3*nabasis*(-(p*(p + 1))/2 + q + nelements*p) + 1 + do l = 1, nbasis3 + ! Get index of where the contributions of atoms i,j,k should be added + z = s + (l - 1)*nabasis + ! Add the contributions for atoms i,j,k + atom_rep(z:z + nabasis - 1) = atom_rep(z:z + nabasis - 1) + angular*radial(l)*rdecay(i, j)*rdecay(i, k) + do t = 1, 3 + ! Add up all gradient contributions wrt atom i + atom_grad(z:z + nabasis - 1, i, t) = atom_grad(z:z + nabasis - 1, i, t) + & + & d_angular*d_angular_d_i(t)*radial(l) + & + & angular*d_radial(l)*d_radial_d_i(t) + & + & angular*radial(l)*(d_ijdecay(t)*rdecay(i, k) + rdecay(i, j)*d_ikdecay(t)) + ! Add up all gradient contributions wrt atom j + atom_grad(z:z + nabasis - 1, j, t) = atom_grad(z:z + nabasis - 1, j, t) + & + & d_angular*d_angular_d_j(t)*radial(l) + & + & angular*d_radial(l)*d_radial_d_j(t) - & + & angular*radial(l)*d_ijdecay(t)*rdecay(i, k) + ! Add up all gradient contributions wrt atom k + atom_grad(z:z + nabasis - 1, k, t) = atom_grad(z:z + nabasis - 1, k, t) + & + & d_angular*d_angular_d_k(t)*radial(l) + & + & angular*d_radial(l)*d_radial_d_k(t) - & + & angular*radial(l)*rdecay(i, j)*d_ikdecay(t) + end do + end do + end do + end do + rep(i, twobody_size + 1:) = atom_rep + grad(i, twobody_size + 1:, :, :) = atom_grad + end do + !$OMP END PARALLEL DO + + deallocate (rdecay) + deallocate (element_types) + deallocate (distance_matrix) + deallocate (inv_distance_matrix) + deallocate (sq_distance_matrix) + deallocate (inv_sq_distance_matrix) + deallocate (atom_rep) + deallocate (atom_grad) + deallocate (a) + deallocate (b) + deallocate (c) + deallocate (radial) + deallocate (angular_base) + deallocate (angular) + deallocate (d_angular) + deallocate (d_angular_d_i) + deallocate (d_angular_d_j) + deallocate (d_angular_d_k) + deallocate (d_radial) + deallocate (d_radial_d_i) + deallocate (d_radial_d_j) + deallocate (d_radial_d_k) + deallocate (d_ijdecay) + deallocate (d_ikdecay) end subroutine fgenerate_acsf_and_gradients - subroutine fgenerate_fchl_acsf(coordinates, nuclear_charges, elements, & - & Rs2, Rs3, Ts, eta2, eta3, zeta, rcut, acut, natoms, natoms_tot, rep_size, & - & two_body_decay, three_body_decay, three_body_weight, rep) - - use acsf_utils, only: decay, calc_angle, calc_cos_angle, true_atom_id - - implicit none - - double precision, intent(in), dimension(:, :) :: coordinates - integer, intent(in), dimension(:) :: nuclear_charges - integer, intent(in), dimension(:) :: elements - double precision, intent(in), dimension(:) :: Rs2 - double precision, intent(in), dimension(:) :: Rs3 - double precision, intent(in), dimension(:) :: Ts - double precision, intent(in) :: eta2 - double precision, intent(in) :: eta3 - double precision, intent(in) :: zeta - double precision, intent(in) :: rcut - double precision, intent(in) :: acut - integer, intent(in) :: natoms, natoms_tot ! natoms_tot includes "virtual" atoms created to satisfy periodic boundary conditions. - integer, intent(in) :: rep_size - double precision, intent(in) :: two_body_decay - double precision, intent(in) :: three_body_decay - double precision, intent(in) :: three_body_weight - - double precision, intent(out), dimension(natoms, rep_size) :: rep - - integer:: two_body_rep_size - integer :: i, j, k, l, n, m, o, p, q, s, z, nelements, nbasis2, nbasis3, nabasis, j_id, j_id_max, true_i - integer, allocatable, dimension(:) :: element_types - double precision :: rij, rik, angle, cos_1, cos_2, cos_3, invcut - ! double precision :: angle_1, angle_2, angle_3 - double precision, allocatable, dimension(:) :: radial, angular, a, b, c - double precision, allocatable, dimension(:, :) :: distance_matrix, rdecay - double precision, allocatable, dimension(:, :) :: saved_radial - integer, allocatable, dimension(:):: saved_j - - double precision :: mu, sigma, ksi3 - - double precision, parameter :: pi = 4.0d0 * atan(1.0d0) - - - if (natoms /= size(nuclear_charges, dim=1)) then - write(*,*) "ERROR: Atom Centered Symmetry Functions creation" - write(*,*) natoms, "coordinates, but", & - & size(nuclear_charges, dim=1), "atom_types!" - stop - endif - - - ! number of element types - nelements = size(elements) - ! Allocate temporary - allocate(element_types(natoms_tot)) - - ! Store element index of every atom - !$OMP PARALLEL DO SCHEDULE(dynamic) - do i = 1, natoms_tot - true_i=true_atom_id(i, natoms) - do j = 1, nelements - if (nuclear_charges(true_i) .eq. elements(j)) then - element_types(i) = j - exit - endif - enddo - enddo - !$OMP END PARALLEL DO - - - ! Get distance matrix - ! Allocate temporary - allocate(distance_matrix(natoms_tot, natoms_tot)) - distance_matrix = 0.0d0 - - - !$OMP PARALLEL DO PRIVATE(rij) SCHEDULE(dynamic) - do i = 1, natoms_tot - do j = i+1, natoms_tot - rij = norm2(coordinates(j,:) - coordinates(i,:)) - distance_matrix(i, j) = rij - distance_matrix(j, i) = rij - enddo - enddo - !$OMP END PARALLEL DO - - ! number of basis functions in the two body term - nbasis2 = size(Rs2) - - ! Inverse of the two body cutoff - invcut = 1.0d0 / rcut - - ! pre-calculate the radial decay in the two body terms - rdecay = decay(distance_matrix, invcut, natoms_tot) - - ! Allocate temporary - allocate(radial(nbasis2)) - allocate(saved_radial(nbasis2, natoms_tot), saved_j(natoms_tot)) - rep=0.0d0 - !$OMP PARALLEL DO PRIVATE(n,m,rij,radial,mu,sigma,saved_radial, saved_j, j_id_max) SCHEDULE(dynamic) - do i = 1, natoms - ! index of the element of atom i - m = element_types(i) - j_id_max=0 - do j = i + 1, natoms_tot - ! distance between atoms i and j - rij = distance_matrix(i,j) - if (rij > rcut) cycle - j_id_max=j_id_max+1 - saved_j(j_id_max)=j - ! index of the element of atom j - n = element_types(j) + & Rs2, Rs3, Ts, eta2, eta3, zeta, rcut, acut, natoms, natoms_tot, rep_size, & + & two_body_decay, three_body_decay, three_body_weight, rep) + + use acsf_utils, only: decay, calc_angle, calc_cos_angle, true_atom_id + + implicit none + + double precision, intent(in), dimension(:, :) :: coordinates + integer, intent(in), dimension(:) :: nuclear_charges + integer, intent(in), dimension(:) :: elements + double precision, intent(in), dimension(:) :: Rs2 + double precision, intent(in), dimension(:) :: Rs3 + double precision, intent(in), dimension(:) :: Ts + double precision, intent(in) :: eta2 + double precision, intent(in) :: eta3 + double precision, intent(in) :: zeta + double precision, intent(in) :: rcut + double precision, intent(in) :: acut + integer, intent(in) :: natoms, natoms_tot ! natoms_tot includes "virtual" atoms created to satisfy periodic boundary conditions. + integer, intent(in) :: rep_size + double precision, intent(in) :: two_body_decay + double precision, intent(in) :: three_body_decay + double precision, intent(in) :: three_body_weight + + double precision, intent(out), dimension(natoms, rep_size) :: rep +! Introduced to make OpenMP parallelization easier. + + integer:: i, j, k, l, n, m, o, p, q, s, z, nelements, nbasis2, nbasis3, nabasis, j_id, max_j_id + integer, allocatable, dimension(:) :: element_types, saved_j + double precision :: rij, rik, angle, cos_1, cos_2, cos_3, invcut + double precision, allocatable, dimension(:) :: radial, angular, a, b, c + double precision, allocatable, dimension(:, :):: saved_radial + double precision, allocatable, dimension(:, :) :: distance_matrix, rdecay + + double precision :: mu, sigma, ksi3 + + double precision, parameter :: pi = 4.0d0*atan(1.0d0) + + if (natoms /= size(nuclear_charges, dim=1)) then + write (*, *) "ERROR: Atom Centered Symmetry Functions creation" + write (*, *) natoms, "coordinates, but", & + & size(nuclear_charges, dim=1), "atom_types!" + stop + end if + +! number of element types + nelements = size(elements) +! Allocate temporary + allocate (element_types(natoms_tot)) + +! Store element index of every atom +!$OMP PARALLEL DO SCHEDULE(dynamic) + do i = 1, natoms_tot + do j = 1, nelements + if (nuclear_charges(true_atom_id(i, natoms)) .eq. elements(j)) then + element_types(i) = j + exit + end if + end do + end do +!$OMP END PARALLEL DO + +! Get distance matrix +! Allocate temporary + allocate (distance_matrix(natoms_tot, natoms_tot)) + distance_matrix = 0.0d0 + +!$OMP PARALLEL DO PRIVATE(rij) SCHEDULE(dynamic) + do i = 1, natoms_tot + do j = i + 1, natoms_tot + rij = norm2(coordinates(j, :) - coordinates(i, :)) + distance_matrix(i, j) = rij + distance_matrix(j, i) = rij + end do + end do +!$OMP END PARALLEL DO + +! number of basis functions in the two body term + nbasis2 = size(Rs2) + +! Inverse of the two body cutoff + invcut = 1.0d0/rcut + +! pre-calculate the radial decay in the two body terms + rdecay = decay(distance_matrix, invcut, natoms_tot) + +! Allocate temporary + allocate (radial(nbasis2), saved_radial(nbasis2, natoms_tot), saved_j(natoms_tot)) + + radial = 0.0d0 +!$OMP PARALLEL DO PRIVATE(n,m,rij,radial,mu,sigma, saved_radial, saved_j) SCHEDULE(dynamic) + do i = 1, natoms +! index of the element of atom i + max_j_id = 0 + do j = i + 1, natoms_tot +! distance between atoms i and j + rij = distance_matrix(i, j) + if (rij <= rcut) then + max_j_id = max_j_id + 1 + saved_j(max_j_id) = j ! two body term of the representation - mu = log(rij / sqrt(1.0d0 + eta2 / rij**2)) - sigma = sqrt(log(1.0d0 + eta2 / rij**2)) + mu = log(rij/sqrt(1.0d0 + eta2/rij**2)) + sigma = sqrt(log(1.0d0 + eta2/rij**2)) radial(:) = 0.0d0 do k = 1, nbasis2 - radial(k) = 1.0d0/(sigma* sqrt(2.0d0*pi) * Rs2(k)) * rdecay(i,j) & - & * exp( - (log(Rs2(k)) - mu)**2 / (2.0d0 * sigma**2) ) / rij**two_body_decay - enddo - saved_radial(:, j_id_max)=radial(:) - enddo - !$OMP CRITICAL - do j_id = 1, j_id_max - j=saved_j(j_id) - n=element_types(j) - rep(i, (n-1)*nbasis2 + 1:n*nbasis2) = rep(i, (n-1)*nbasis2 + 1:n*nbasis2) + saved_radial(:, j_id) - if (j<=natoms) rep(j, (m-1)*nbasis2 + 1:m*nbasis2) = rep(j, (m-1)*nbasis2 + 1:m*nbasis2) + saved_radial(:, j_id) - enddo - !$OMP END CRITICAL - enddo - !$OMP END PARALLEL DO - - deallocate(radial) - deallocate(saved_radial, saved_j) - - ! number of radial basis functions in the three body term - nbasis3 = size(Rs3) - ! number of radial basis functions in the three body term - nabasis = size(Ts) - - ! Inverse of the three body cutoff - invcut = 1.0d0 / acut - ! pre-calculate the radial decay in the three body terms - rdecay = decay(distance_matrix, invcut, natoms_tot) - - ! Allocate temporary - allocate(a(3)) - allocate(b(3)) - allocate(c(3)) - allocate(radial(nbasis3)) - allocate(angular(nabasis)) - two_body_rep_size=nbasis2*nelements - - ! This could probably be done more efficiently if it's a bottleneck - ! Also the order is a bit wobbly compared to the tensorflow implementation - !$OMP PARALLEL DO PRIVATE(rij, n, rik, m, a, b, c, angle, radial, angular, & - !$OMP cos_1, cos_2, cos_3, mu, sigma, o, ksi3, & - !$OMP p, q, s, z, l) SCHEDULE(dynamic) - do i = 1, natoms - do j = 1, natoms_tot - 1 - if (i .eq. j) cycle - ! distance between atoms i and j - rij = distance_matrix(i,j) - if (rij > acut) cycle - ! index of the element of atom j - n = element_types(j) - do k = j + 1, natoms_tot - if (i .eq. k) cycle - if (j .eq. k) cycle - ! distance between atoms i and k - rik = distance_matrix(i,k) - if (rik > acut) cycle - ! index of the element of atom k - m = element_types(k) - ! coordinates of atoms j, i, k - a = coordinates(j,:) - b = coordinates(i,:) - c = coordinates(k,:) - ! angle between atoms i, j and k centered on i - angle = calc_angle(a,b,c) - cos_1 = calc_cos_angle(a,b,c) - cos_2 = calc_cos_angle(a,c,b) - cos_3 = calc_cos_angle(b,a,c) - - ! The radial part of the three body terms including decay - radial = exp(-eta3*(0.5d0 * (rij+rik) - Rs3)**2) * rdecay(i,j) * rdecay(i,k) - - ksi3 = (1.0d0 + 3.0d0 * cos_1 * cos_2 * cos_3) & - & / (distance_matrix(i,k) * distance_matrix(i,j) * distance_matrix(j,k) & - & )**three_body_decay * three_body_weight - - angular = 0.0d0 - do l = 1, nabasis/2 - - o = l*2-1 - angular(2*l-1) = angular(2*l-1) + 2*cos(o * angle) & - & * exp(-(zeta * o)**2 /2) - - angular(2*l) = angular(2*l) + 2*sin(o * angle) & - & * exp(-(zeta * o)**2 /2) - - enddo - - ! The lowest of the element indices for atoms j and k - p = min(n,m) - 1 - ! The highest of the element indices for atoms j and k - q = max(n,m) - 1 - ! calculate the indices that the three body terms should be added to - s = two_body_rep_size + nbasis3 * nabasis * (-(p * (p + 1))/2 + q + nelements * p) + 1 - - do l = 1, nbasis3 - ! calculate the indices that the three body terms should be added to - z = s + (l-1) * nabasis - ! Add the contributions from atoms i,j and k - rep(i, z:z + nabasis - 1) = rep(i, z:z + nabasis - 1) + angular * radial(l) * ksi3 - enddo - enddo - enddo - enddo - !$OMP END PARALLEL DO - - deallocate(element_types) - deallocate(rdecay) - deallocate(distance_matrix) - deallocate(a) - deallocate(b) - deallocate(c) - deallocate(radial) - deallocate(angular) + radial(k) = 1.0d0/(sigma*sqrt(2.0d0*pi)*Rs2(k))*rdecay(i, j) & + & *exp(-(log(Rs2(k)) - mu)**2/(2.0d0*sigma**2))/rij**two_body_decay + end do + saved_radial(:, max_j_id) = radial(:) + end if + end do + m = element_types(i) + !$OMP CRITICAL + do j_id = 1, max_j_id + j = saved_j(j_id) + ! index of the element of atom j + n = element_types(j) + rep(i, (n - 1)*nbasis2 + 1:n*nbasis2) = rep(i, (n - 1)*nbasis2 + 1:n*nbasis2) + saved_radial(:, j_id) + if (j <= natoms) then + rep(j, (m - 1)*nbasis2 + 1:m*nbasis2) = rep(j, (m - 1)*nbasis2 + 1:m*nbasis2) + saved_radial(:, j_id) + end if + end do + !$OMP END CRITICAL + end do +!$OMP END PARALLEL DO + deallocate (radial, saved_radial, saved_j) + +! number of radial basis functions in the three body term + nbasis3 = size(Rs3) +! number of radial basis functions in the three body term + nabasis = size(Ts) + +! Inverse of the three body cutoff + invcut = 1.0d0/acut +! pre-calculate the radial decay in the three body terms + rdecay = decay(distance_matrix, invcut, natoms_tot) + +! Allocate temporary + allocate (a(3)) + allocate (b(3)) + allocate (c(3)) + allocate (radial(nbasis3)) + allocate (angular(nabasis)) + +!$OMP PARALLEL DO PRIVATE(rij, n, rik, m, a, b, c, angle, cos_1, cos_2, cos_3,& +!$OMP radial, ksi3, angular, o, p, q, s, z) SCHEDULE(dynamic) + do i = 1, natoms + do j = 1, natoms_tot - 1 + if (i .eq. j) cycle +! distance between atoms i and j + rij = distance_matrix(i, j) + if (rij > acut) cycle +! index of the element of atom j + n = element_types(j) + do k = j + 1, natoms_tot + if (i .eq. k) cycle + if (j .eq. k) cycle +! distance between atoms i and k + rik = distance_matrix(i, k) + if (rik > acut) cycle +! index of the element of atom k + m = element_types(k) +! coordinates of atoms j, i, k + a = coordinates(j, :) + b = coordinates(i, :) + c = coordinates(k, :) +! angle between atoms i, j and k centered on i + angle = calc_angle(a, b, c) + cos_1 = calc_cos_angle(a, b, c) + cos_2 = calc_cos_angle(a, c, b) + cos_3 = calc_cos_angle(b, a, c) + +! The radial part of the three body terms including decay + radial = exp(-eta3*(0.5d0*(rij + rik) - Rs3)**2)*rdecay(i, j)*rdecay(i, k) + + ksi3 = (1.0d0 + 3.0d0*cos_1*cos_2*cos_3) & + & /(distance_matrix(i, k)*distance_matrix(i, j)*distance_matrix(j, k) & + & )**three_body_decay*three_body_weight + + angular = 0.0d0 + do l = 1, nabasis/2 + + o = l*2 - 1 + angular(2*l - 1) = angular(2*l - 1) + 2*cos(o*angle) & + & *exp(-(zeta*o)**2/2) + + angular(2*l) = angular(2*l) + 2*sin(o*angle) & + & *exp(-(zeta*o)**2/2) + + end do + +! The lowest of the element indices for atoms j and k + p = min(n, m) - 1 +! The highest of the element indices for atoms j and k + q = max(n, m) - 1 +! calculate the indices that the three body terms should be added to + s = nelements*nbasis2 + nbasis3*nabasis*(-(p*(p + 1))/2 + q + nelements*p) + 1 + + do l = 1, nbasis3 +! calculate the indices that the three body terms should be added to + z = s + (l - 1)*nabasis +! Add the contributions from atoms i,j and k + rep(i, z:z + nabasis - 1) = rep(i, z:z + nabasis - 1) + angular*radial(l)*ksi3 + end do + end do + end do + end do +!$OMP END PARALLEL DO + + deallocate (element_types) + deallocate (rdecay) + deallocate (distance_matrix) + deallocate (a) + deallocate (b) + deallocate (c) + deallocate (radial) + deallocate (angular) end subroutine fgenerate_fchl_acsf subroutine fgenerate_fchl_acsf_and_gradients(coordinates, nuclear_charges, elements, & - & Rs2, Rs3, Ts, eta2, eta3, zeta, rcut, acut, natoms, natoms_tot, rep_size, & - & two_body_decay, three_body_decay, three_body_weight, rep, grad) - - use acsf_utils, only: decay, calc_angle, calc_cos_angle, true_atom_id - - implicit none - - double precision, intent(in), dimension(:, :) :: coordinates - integer, intent(in), dimension(:) :: nuclear_charges - integer, intent(in), dimension(:) :: elements - double precision, intent(in), dimension(:) :: Rs2 - double precision, intent(in), dimension(:) :: Rs3 - double precision, intent(in), dimension(:) :: Ts - double precision, intent(in) :: eta2 - double precision, intent(in) :: eta3 - double precision, intent(in) :: zeta - double precision, intent(in) :: rcut - double precision, intent(in) :: acut - - double precision, intent(in) :: two_body_decay - double precision, intent(in) :: three_body_decay - double precision, intent(in) :: three_body_weight - - double precision :: mu, sigma, dx, exp_s2, scaling, dscal, ddecay - double precision :: cos_i, cos_j, cos_k - double precision, allocatable, dimension(:) :: exp_ln - double precision, allocatable, dimension(:) :: log_Rs2 - - integer, intent(in) :: natoms, natoms_tot ! natoms_tot includes "virtual" atoms created to satisfy periodic boundary conditions. - integer, intent(in) :: rep_size - double precision, intent(out), dimension(natoms, rep_size) :: rep - double precision, intent(out), dimension(natoms, rep_size, natoms, 3) :: grad + & Rs2, Rs3, Ts, eta2, eta3, zeta, rcut, acut, natoms, natoms_tot, rep_size, & + & two_body_decay, three_body_decay, three_body_weight, rep, grad) + + use acsf_utils, only: decay, calc_angle, calc_cos_angle, true_atom_id + + implicit none + + double precision, intent(in), dimension(:, :) :: coordinates + integer, intent(in), dimension(:) :: nuclear_charges + integer, intent(in), dimension(:) :: elements + double precision, intent(in), dimension(:) :: Rs2 + double precision, intent(in), dimension(:) :: Rs3 + double precision, intent(in), dimension(:) :: Ts + double precision, intent(in) :: eta2 + double precision, intent(in) :: eta3 + double precision, intent(in) :: zeta + double precision, intent(in) :: rcut + double precision, intent(in) :: acut + + double precision, intent(in) :: two_body_decay + double precision, intent(in) :: three_body_decay + double precision, intent(in) :: three_body_weight + + double precision :: mu, sigma, dx, exp_s2, scaling, dscal, ddecay + double precision :: cos_i, cos_j, cos_k + double precision, allocatable, dimension(:) :: exp_ln + double precision, allocatable, dimension(:) :: log_Rs2 + + integer, intent(in) :: natoms, natoms_tot ! natoms_tot includes "virtual" atoms created to satisfy periodic boundary conditions. + integer, intent(in) :: rep_size + double precision, intent(out), dimension(natoms, rep_size) :: rep + double precision, intent(out), dimension(natoms, rep_size, natoms, 3) :: grad ! Introduced to make OpenMP parallelization easier. - double precision, dimension(:, :, :), allocatable:: add_rep - double precision, dimension(:, :, :, :), allocatable :: add_grad - - integer :: i, j, k, l, m, n, p, q, s, t, z, nelements, nbasis2, nbasis3, nabasis, twobody_size, true_j, true_k - integer, allocatable, dimension(:) :: element_types - double precision :: rij, rik, angle, dot, rij2, rik2, invrij, invrik, invrij2, invrik2, invcut - double precision, allocatable, dimension(:) :: radial_base, radial, angular, a, b, c, atom_rep - double precision, allocatable, dimension(:) :: angular_base, d_radial, d_radial_d_i - double precision, allocatable, dimension(:) :: d_radial_d_j, d_radial_d_k, d_ijdecay, d_ikdecay - double precision, allocatable, dimension(:) :: d_angular, part, radial_part - double precision, allocatable, dimension(:) :: d_angular_d_i, d_angular_d_j, d_angular_d_k - double precision, allocatable, dimension(:, :) :: distance_matrix, rdecay, sq_distance_matrix - double precision, allocatable, dimension(:, :) :: inv_distance_matrix, inv_sq_distance_matrix - double precision, allocatable, dimension(:, :, :) :: atom_grad - - double precision :: atm, atm_i, atm_j, atm_k - double precision :: invrjk, invr_atm, vi, vj, vk - double precision, allocatable, dimension(:) :: d_atm_i, d_atm_j, d_atm_k - double precision, allocatable, dimension(:) :: d_atm_ii, d_atm_ji, d_atm_ki - double precision, allocatable, dimension(:) :: d_atm_ij, d_atm_jj, d_atm_kj - double precision, allocatable, dimension(:) :: d_atm_ik, d_atm_jk, d_atm_kk - double precision, allocatable, dimension(:) :: d_atm_i2, d_atm_j2, d_atm_k2 - double precision, allocatable, dimension(:) :: d_atm_i3, d_atm_j3, d_atm_k3 - double precision, allocatable, dimension(:) :: d_atm_extra_i, d_atm_extra_j, d_atm_extra_k - - double precision, parameter :: pi = 4.0d0 * atan(1.0d0) - - if (natoms /= size(nuclear_charges, dim=1)) then - write(*,*) "ERROR: Atom Centered Symmetry Functions creation" - write(*,*) natoms, "coordinates, but", & - & size(nuclear_charges, dim=1), "atom_types!" - stop - endif - - ! Number of unique elements - nelements = size(elements) - ! Allocate temporary - allocate(element_types(natoms)) - - ! Store element index of every atom - !$OMP PARALLEL DO SCHEDULE(dynamic) - do i = 1, natoms - do j = 1, nelements - if (nuclear_charges(i) .eq. elements(j)) then - element_types(i) = j - exit - endif - enddo - enddo - !$OMP END PARALLEL DO - - - - ! Get distance matrix - ! Allocate temporary - allocate(distance_matrix(natoms_tot, natoms_tot)) - allocate(sq_distance_matrix(natoms_tot, natoms_tot)) - allocate(inv_distance_matrix(natoms_tot, natoms_tot)) - allocate(inv_sq_distance_matrix(natoms_tot, natoms_tot)) - distance_matrix = 0.0d0 - sq_distance_matrix = 0.0d0 - inv_distance_matrix = 0.0d0 - inv_sq_distance_matrix = 0.0d0 - - - !$OMP PARALLEL DO PRIVATE(rij,rij2,invrij,invrij2) SCHEDULE(dynamic) - do i = 1, natoms_tot - do j = i+1, natoms_tot - rij = norm2(coordinates(j,:) - coordinates(i,:)) - distance_matrix(i, j) = rij - distance_matrix(j, i) = rij - rij2 = rij * rij - sq_distance_matrix(i, j) = rij2 - sq_distance_matrix(j, i) = rij2 - invrij = 1.0d0 / rij - inv_distance_matrix(i, j) = invrij - inv_distance_matrix(j, i) = invrij - invrij2 = invrij * invrij - inv_sq_distance_matrix(i, j) = invrij2 - inv_sq_distance_matrix(j, i) = invrij2 - enddo - enddo - !$OMP END PARALLEL DO - - - ! Number of two body basis functions - nbasis2 = size(Rs2) - - ! Inverse of the two body cutoff distance - invcut = 1.0d0 / rcut - ! Pre-calculate the two body decay - rdecay = decay(distance_matrix, invcut, natoms_tot) - - ! Allocate temporary - allocate(radial_base(nbasis2)) - allocate(radial(nbasis2)) - allocate(radial_part(nbasis2)) - allocate(part(nbasis2)) - allocate(exp_ln(nbasis2)) - allocate(log_Rs2(nbasis2)) - - grad =0.0d0 - rep = 0.0d0 - - allocate(add_rep(nbasis2, natoms, natoms), add_grad(nbasis2, 3, natoms, natoms)) - add_rep=0.0d0 - add_grad=0.0d0 - log_Rs2(:) = log(Rs2(:)) - - !$OMP PARALLEL DO PRIVATE(m, n, rij, invrij, mu, sigma, exp_s2, exp_ln,& - !$OMP scaling, radial_base, radial, dx, part, dscal, ddecay) SCHEDULE(dynamic) - do i = 1, natoms - ! The element index of atom i - m = element_types(i) - do j = i + 1, natoms_tot - rij = distance_matrix(i,j) - if (rij > rcut) continue - true_j=true_atom_id(i, natoms) - if (true_j < i) continue - ! The element index of atom j - n = element_types(true_j) - ! Distance between atoms i and j - invrij = inv_distance_matrix(i,j) - mu = log(rij / sqrt(1.0d0 + eta2 * inv_sq_distance_matrix(i, j))) - sigma = sqrt(log(1.0d0 + eta2 * inv_sq_distance_matrix(i, j))) + double precision, dimension(:, :, :), allocatable:: add_rep + double precision, dimension(:, :, :, :), allocatable :: add_grad + + integer :: i, j, k, l, m, n, p, q, s, t, z, nelements, nbasis2, nbasis3, nabasis, twobody_size, true_j, true_k + integer, allocatable, dimension(:) :: element_types + double precision :: rij, rik, angle, dot, rij2, rik2, invrij, invrik, invrij2, invrik2, invcut + double precision, allocatable, dimension(:) :: radial_base, radial, angular, a, b, c, atom_rep + double precision, allocatable, dimension(:) :: angular_base, d_radial, d_radial_d_i + double precision, allocatable, dimension(:) :: d_radial_d_j, d_radial_d_k, d_ijdecay, d_ikdecay + double precision, allocatable, dimension(:) :: d_angular, part, radial_part + double precision, allocatable, dimension(:) :: d_angular_d_i, d_angular_d_j, d_angular_d_k + double precision, allocatable, dimension(:, :) :: distance_matrix, rdecay, sq_distance_matrix + double precision, allocatable, dimension(:, :) :: inv_distance_matrix, inv_sq_distance_matrix + double precision, allocatable, dimension(:, :, :) :: atom_grad + + double precision :: atm, atm_i, atm_j, atm_k + double precision :: invrjk, invr_atm, vi, vj, vk + double precision, allocatable, dimension(:) :: d_atm_i, d_atm_j, d_atm_k + double precision, allocatable, dimension(:) :: d_atm_ii, d_atm_ji, d_atm_ki + double precision, allocatable, dimension(:) :: d_atm_ij, d_atm_jj, d_atm_kj + double precision, allocatable, dimension(:) :: d_atm_ik, d_atm_jk, d_atm_kk + double precision, allocatable, dimension(:) :: d_atm_i2, d_atm_j2, d_atm_k2 + double precision, allocatable, dimension(:) :: d_atm_i3, d_atm_j3, d_atm_k3 + double precision, allocatable, dimension(:) :: d_atm_extra_i, d_atm_extra_j, d_atm_extra_k + + double precision, parameter :: pi = 4.0d0*atan(1.0d0) + + if (natoms /= size(nuclear_charges, dim=1)) then + write (*, *) "ERROR: Atom Centered Symmetry Functions creation" + write (*, *) natoms, "coordinates, but", & + & size(nuclear_charges, dim=1), "atom_types!" + stop + end if + + ! Number of unique elements + nelements = size(elements) + ! Allocate temporary + allocate (element_types(natoms)) + +! Store element index of every atom +!$OMP PARALLEL DO SCHEDULE(dynamic) + do i = 1, natoms + do j = 1, nelements + if (nuclear_charges(i) .eq. elements(j)) then + element_types(i) = j + exit + end if + end do + end do +!$OMP END PARALLEL DO + + ! Get distance matrix + ! Allocate temporary + allocate (distance_matrix(natoms_tot, natoms_tot)) + allocate (sq_distance_matrix(natoms_tot, natoms_tot)) + allocate (inv_distance_matrix(natoms_tot, natoms_tot)) + allocate (inv_sq_distance_matrix(natoms_tot, natoms_tot)) + distance_matrix = 0.0d0 + sq_distance_matrix = 0.0d0 + inv_distance_matrix = 0.0d0 + inv_sq_distance_matrix = 0.0d0 + +!$OMP PARALLEL DO PRIVATE(rij,rij2,invrij,invrij2) SCHEDULE(dynamic) + do i = 1, natoms_tot + do j = i + 1, natoms_tot + rij = norm2(coordinates(j, :) - coordinates(i, :)) + distance_matrix(i, j) = rij + distance_matrix(j, i) = rij + rij2 = rij*rij + sq_distance_matrix(i, j) = rij2 + sq_distance_matrix(j, i) = rij2 + invrij = 1.0d0/rij + inv_distance_matrix(i, j) = invrij + inv_distance_matrix(j, i) = invrij + invrij2 = invrij*invrij + inv_sq_distance_matrix(i, j) = invrij2 + inv_sq_distance_matrix(j, i) = invrij2 + end do + end do +!$OMP END PARALLEL DO + + ! Number of two body basis functions + nbasis2 = size(Rs2) + + ! Inverse of the two body cutoff distance + invcut = 1.0d0/rcut + ! Pre-calculate the two body decay + rdecay = decay(distance_matrix, invcut, natoms_tot) + + ! Allocate temporary + allocate (radial_base(nbasis2)) + allocate (radial(nbasis2)) + allocate (radial_part(nbasis2)) + allocate (part(nbasis2)) + allocate (exp_ln(nbasis2)) + allocate (log_Rs2(nbasis2)) + + grad = 0.0d0 + rep = 0.0d0 + + allocate (add_rep(nbasis2, natoms, natoms), add_grad(nbasis2, 3, natoms, natoms)) + add_rep = 0.0d0 + add_grad = 0.0d0 + log_Rs2(:) = log(Rs2(:)) + +!$OMP PARALLEL DO PRIVATE(m, n, rij, invrij, mu, sigma, exp_s2, exp_ln,& +!$OMP scaling, radial_base, radial, dx, part, dscal, ddecay) SCHEDULE(dynamic) + do i = 1, natoms + ! The element index of atom i + m = element_types(i) + do j = i + 1, natoms_tot + ! The element index of atom j + true_j = true_atom_id(j, natoms) + n = element_types(true_j) + ! Distance between atoms i and j + rij = distance_matrix(i, j) + if (rij <= rcut) then + invrij = inv_distance_matrix(i, j) + + mu = log(rij/sqrt(1.0d0 + eta2*inv_sq_distance_matrix(i, j))) + sigma = sqrt(log(1.0d0 + eta2*inv_sq_distance_matrix(i, j))) exp_s2 = exp(sigma**2) - exp_ln = exp(-(log_Rs2(:) - mu)**2 / sigma**2 * 0.5d0) * sqrt(2.0d0) + exp_ln = exp(-(log_Rs2(:) - mu)**2/sigma**2*0.5d0)*sqrt(2.0d0) - scaling = 1.0d0 / rij**two_body_decay + scaling = 1.0d0/rij**two_body_decay + radial_base(:) = 1.0d0/(sigma*sqrt(2.0d0*pi)*Rs2(:))*exp(-(log_Rs2(:) - mu)**2/(2.0d0*sigma**2)) - radial_base(:) = 1.0d0/(sigma* sqrt(2.0d0*pi) * Rs2(:)) * exp(-(log_Rs2(:) - mu)**2 / (2.0d0 * sigma**2)) + radial(:) = radial_base(:)*scaling*rdecay(i, j) - radial(:) = radial_base(:) * scaling * rdecay(i,j) + rep(i, (n - 1)*nbasis2 + 1:n*nbasis2) = rep(i, (n - 1)*nbasis2 + 1:n*nbasis2) + radial + if (j <= natoms) add_rep(:, i, j) = radial - rep(i, (n-1)*nbasis2 + 1:n*nbasis2) = rep(i, (n-1)*nbasis2 + 1:n*nbasis2) + radial - add_rep(:, i, true_j)=add_rep(:, i, true_j)+radial do k = 1, 3 - dx = -(coordinates(i,k) - coordinates(j,k)) - part(:) = ((log_Rs2(:) - mu) * (-dx *(rij**2 * exp_s2 + eta2) / (rij * sqrt(exp_s2))**3) & - &* sqrt(exp_s2) / (sigma**2 * rij) + (log_Rs2(:) - mu) ** 2 * eta2 * dx / & - &(sigma**4 * rij**4 * exp_s2)) * exp_ln / (Rs2(:) * sigma * sqrt(pi) * 2) & - &- exp_ln * eta2 * dx / (Rs2(:) * sigma**3 *sqrt(pi) * rij**4 * exp_s2 * 2.0d0) - - dscal = two_body_decay * dx / rij**(two_body_decay+2.0d0) - ddecay = dx * 0.5d0 * pi * sin(pi*rij * invcut) * invcut * invrij - - part(:) = part(:) * scaling * rdecay(i,j) + radial_base(:) * dscal * rdecay(i,j) & - & + radial_base(:) * scaling * ddecay - - ! The gradients wrt coordinates - grad(i, (n-1)*nbasis2 + 1:n*nbasis2, i, k) = grad(i, (n-1)*nbasis2 + 1:n*nbasis2, i, k) + part - grad(i, (n-1)*nbasis2 + 1:n*nbasis2, true_j, k) = grad(i, (n-1)*nbasis2 + 1:n*nbasis2, true_j, k) - part - grad(true_j, (m-1)*nbasis2 + 1:m*nbasis2, i, k) = grad(true_j, (m-1)*nbasis2 + 1:m*nbasis2, i, k) + part - add_grad(:, k, i, true_j)=add_grad(:, k, i, true_j)+part - enddo - enddo - enddo - !$OMP END PARALLEL DO - - !$OMP PARALLEL DO PRIVATE(m) SCHEDULE(dynamic) - do j = 2, natoms - do i = 1, j-1 - if (distance_matrix(i,j)>rcut) cycle - m = element_types(i) - rep(j, (m-1)*nbasis2 + 1:m*nbasis2) = rep(j, (m-1)*nbasis2 + 1:m*nbasis2) + add_rep(:, i, j) - do k=1, 3 - grad(j, (m-1)*nbasis2 + 1:m*nbasis2, j, k) = grad(j, (m-1)*nbasis2 + 1:m*nbasis2, j, k) - add_grad(:, k, i, j) - enddo - enddo - enddo - !$OMP END PARALLEL DO - - deallocate(add_rep, add_grad) - - deallocate(radial_base) - deallocate(radial) - deallocate(radial_part) - deallocate(part) - - - ! Number of radial basis functions in the three body term - nbasis3 = size(Rs3) - ! Number of angular basis functions in the three body term - nabasis = size(Ts) - ! Size of two body terms - twobody_size = nelements * nbasis2 - - ! Inverse of the three body cutoff distance - invcut = 1.0d0 / acut - ! Pre-calculate the three body decay - rdecay = decay(distance_matrix, invcut, natoms_tot) - - ! Allocate temporary - allocate(atom_rep(rep_size - twobody_size)) - allocate(atom_grad(rep_size - twobody_size, natoms, 3)) - allocate(a(3)) - allocate(b(3)) - allocate(c(3)) - allocate(radial(nbasis3)) - allocate(radial_base(nbasis3)) - allocate(angular_base(nabasis)) - allocate(angular(nabasis)) - allocate(d_angular(nabasis)) - allocate(d_angular_d_i(3)) - allocate(d_angular_d_j(3)) - allocate(d_angular_d_k(3)) - allocate(d_radial(nbasis3)) - allocate(d_radial_d_i(3)) - allocate(d_radial_d_j(3)) - allocate(d_radial_d_k(3)) - allocate(d_ijdecay(3)) - allocate(d_ikdecay(3)) - - allocate(d_atm_i(3)) - allocate(d_atm_j(3)) - allocate(d_atm_k(3)) - allocate(d_atm_ii(3)) - allocate(d_atm_ij(3)) - allocate(d_atm_ik(3)) - allocate(d_atm_ji(3)) - allocate(d_atm_jj(3)) - allocate(d_atm_jk(3)) - allocate(d_atm_ki(3)) - allocate(d_atm_kj(3)) - allocate(d_atm_kk(3)) - allocate(d_atm_i2(3)) - allocate(d_atm_j2(3)) - allocate(d_atm_k2(3)) - allocate(d_atm_i3(3)) - allocate(d_atm_j3(3)) - allocate(d_atm_k3(3)) - allocate(d_atm_extra_i(3)) - allocate(d_atm_extra_j(3)) - allocate(d_atm_extra_k(3)) - - !$OMP PARALLEL DO PRIVATE(atom_rep, atom_grad, rij, n, rij2, invrij, invrij2,& - !$OMP rik, m, rik2, invrik, invrjk, invrik2, a, b, c, angle, cos_i, cos_k,& - !$OMP cos_j, radial_base, radial, p, q, dot, angular, d_angular, d_angular_d_j,& - !$OMP d_angular_d_k, d_angular_d_i, d_radial, d_radial_d_j, d_radial_d_k,& - !$OMP d_radial_d_i, d_ijdecay, d_ikdecay, invr_atm, atm, atm_i, atm_j, atm_k, vi,& - !$OMP vj, vk, d_atm_ii, d_atm_ij, d_atm_ik, d_atm_ji, d_atm_jj, d_atm_jk, d_atm_ki,& - !$OMP d_atm_kj, d_atm_kk, d_atm_extra_i, d_atm_extra_j, d_atm_extra_k, s, z) SCHEDULE(dynamic) - do i = 1, natoms - atom_rep = 0.0d0 - atom_grad = 0.0d0 - do j = 1, natoms_tot - 1 - if (i .eq. j) cycle - ! distance between atom i and j - rij = distance_matrix(i,j) - if (rij > acut) cycle - ! index of the element of atom j - true_j=true_atom_id(j, natoms) - n = element_types(true_j) - ! squared distance between atom i and j - rij2 = sq_distance_matrix(i,j) - ! inverse distance between atom i and j - invrij = inv_distance_matrix(i,j) - ! inverse squared distance between atom i and j - invrij2 = inv_sq_distance_matrix(i,j) - do k = j + 1, natoms_tot - if (i .eq. k) cycle - ! distance between atom i and k - rik = distance_matrix(i,k) - if (rik > acut) cycle - ! index of the element of atom k - true_k=true_atom_id(k, natoms) - m = element_types(true_k) - ! squared distance between atom i and k - rik2 = sq_distance_matrix(i,k) - ! inverse distance between atom i and k - invrik = inv_distance_matrix(i,k) - ! inverse distance between atom j and k - invrjk = inv_distance_matrix(j,k) - ! inverse squared distance between atom i and k - invrik2 = inv_sq_distance_matrix(i,k) - ! coordinates of atoms j, i, k - a = coordinates(j,:) - b = coordinates(i,:) - c = coordinates(k,:) - ! angle between atom i, j and k, centered on i - angle = calc_angle(a,b,c) - cos_i = calc_cos_angle(a,b,c) - cos_k = calc_cos_angle(a,c,b) - cos_j = calc_cos_angle(b,a,c) - - ! part of the radial part of the 3body terms - radial_base(:) = exp(-eta3*(0.5d0 * (rij+rik) - Rs3(:))**2) - radial(:) = radial_base(:) ! * scaling - - p = min(n,m) - 1 - ! the highest index of the elements of j,k - q = max(n,m) - 1 - ! Dot product between the vectors connecting atom i,j and i,k - dot = dot_product(a-b,c-b) - - angular(1) = exp(-(zeta**2)*0.5d0) * 2 * cos(angle) - angular(2) = exp(-(zeta**2)*0.5d0) * 2 * sin(angle) - - d_angular(1) = exp(-(zeta**2)*0.5d0) * 2 * sin(angle) / sqrt(max(1d-10, rij2 * rik2 - dot**2)) - d_angular(2) = -exp(-(zeta**2)*0.5d0) * 2 * cos(angle) / sqrt(max(1d-10, rij2 * rik2 - dot**2)) - - ! Part of the derivative of the angular basis functions wrt atom j (dim(3)) - d_angular_d_j = c - b + dot * ((b - a) * invrij2) - ! Part of the derivative of the angular basis functions wrt atom k (dim(3)) - d_angular_d_k = a - b + dot * ((b - c) * invrik2) - ! Part of the derivative of the angular basis functions wrt atom i (dim(3)) - d_angular_d_i = - (d_angular_d_j + d_angular_d_k) - - ! Part of the derivative of the radial basis functions wrt coordinates (dim(nbasis3)) - ! including decay - d_radial = radial * eta3 * (0.5d0 * (rij+rik) - Rs3) ! * rdecay(i,j) * rdecay(i,k) - ! Part of the derivative of the radial basis functions wrt atom j (dim(3)) - d_radial_d_j = (b - a) * invrij - ! Part of the derivative of the radial basis functions wrt atom k (dim(3)) - d_radial_d_k = (b - c) * invrik - ! Part of the derivative of the radial basis functions wrt atom i (dim(3)) - d_radial_d_i = - (d_radial_d_j + d_radial_d_k) - - ! Part of the derivative of the i,j decay functions wrt coordinates (dim(3)) - d_ijdecay = - pi * (b - a) * sin(pi * rij * invcut) * 0.5d0 * invrij * invcut - ! Part of the derivative of the i,k decay functions wrt coordinates (dim(3)) - d_ikdecay = - pi * (b - c) * sin(pi * rik * invcut) * 0.5d0 * invrik * invcut - - invr_atm = (invrij * invrjk *invrik)**three_body_decay - - ! Axilrod-Teller-Muto term - atm = (1.0d0 + 3.0d0 * cos_i * cos_j * cos_k) * invr_atm * three_body_weight - - atm_i = (3.0d0 * cos_j * cos_k) * invr_atm * invrij * invrik - atm_j = (3.0d0 * cos_k * cos_i) * invr_atm * invrij * invrjk - atm_k = (3.0d0 * cos_i * cos_j) * invr_atm * invrjk * invrik - - vi = dot_product(a-b,c-b) - vj = dot_product(c-a,b-a) - vk = dot_product(b-c,a-c) - - d_atm_ii(:) = 2 * b - a - c - vi * ((b-a)*invrij**2 + (b-c)*invrik**2) - d_atm_ij(:) = c - a - vj * (b-a)*invrij**2 - d_atm_ik(:) = a - c - vk * (b-c)*invrik**2 - - d_atm_ji(:) = c - b - vi * (a-b)*invrij**2 - d_atm_jj(:) = 2 * a - b - c - vj * ((a-b)*invrij**2 + (a-c)*invrjk**2) - d_atm_jk(:) = b - c - vk * (a-c)*invrjk**2 - - d_atm_ki(:) = a - b - vi * (c-b)*invrik**2 - d_atm_kj(:) = b - a - vj * (c-a)*invrjk**2 - d_atm_kk(:) = 2 * c - a - b - vk * ((c-a)*invrjk**2 + (c-b)*invrik**2) - - d_atm_extra_i(:) = ((a-b)*invrij**2 + (c-b)*invrik**2) * atm * three_body_decay / three_body_weight - d_atm_extra_j(:) = ((b-a)*invrij**2 + (c-a)*invrjk**2) * atm * three_body_decay / three_body_weight - d_atm_extra_k(:) = ((a-c)*invrjk**2 + (b-c)*invrik**2) * atm * three_body_decay / three_body_weight - - ! Get index of where the contributions of atoms i,j,k should be added - s = nbasis3 * nabasis * (-(p * (p + 1))/2 + q + nelements * p) + 1 - - - do l = 1, nbasis3 - - ! Get index of where the contributions of atoms i,j,k should be added - z = s + (l-1) * nabasis - - ! Add the contributions for atoms i,j,k - atom_rep(z:z + nabasis - 1) = atom_rep(z:z + nabasis - 1) & - & + angular * radial(l) * atm * rdecay(i,j) * rdecay(i,k) - - do t = 1, 3 - - ! Add up all gradient contributions wrt atom i - atom_grad(z:z + nabasis - 1, i, t) = atom_grad(z:z + nabasis - 1, i, t) + & - & d_angular * d_angular_d_i(t) * radial(l) * atm * rdecay(i,j) * rdecay(i,k) + & - & angular * d_radial(l) * d_radial_d_i(t) * atm * rdecay(i,j) * rdecay(i,k) + & - & angular * radial(l) * (atm_i * d_atm_ii(t) + atm_j * d_atm_ij(t) & - & + atm_k * d_atm_ik(t) + d_atm_extra_i(t)) * three_body_weight * rdecay(i,j) * rdecay(i,k) + & - & angular * radial(l) * (d_ijdecay(t) * rdecay(i,k) + rdecay(i,j) * d_ikdecay(t)) * atm - ! Add up all gradient contributions wrt atom j - atom_grad(z:z + nabasis - 1, true_j, t) = atom_grad(z:z + nabasis - 1, true_j, t) + & - & d_angular * d_angular_d_j(t) * radial(l) * atm * rdecay(i,j) * rdecay(i,k) + & - & angular * d_radial(l) * d_radial_d_j(t) * atm * rdecay(i,j) * rdecay(i,k) + & - & angular * radial(l) * (atm_i * d_atm_ji(t) + atm_j * d_atm_jj(t) & - & + atm_k * d_atm_jk(t) + d_atm_extra_j(t)) * three_body_weight * rdecay(i,j) * rdecay(i,k) - & - & angular * radial(l) * d_ijdecay(t) * rdecay(i,k) * atm - - ! Add up all gradient contributions wrt atom k - atom_grad(z:z + nabasis - 1, true_k, t) = atom_grad(z:z + nabasis - 1, true_k, t) + & - & d_angular * d_angular_d_k(t) * radial(l) * atm * rdecay(i,j) * rdecay(i,k) + & - & angular * d_radial(l) * d_radial_d_k(t) * atm * rdecay(i,j) * rdecay(i,k) + & - & angular * radial(l) * (atm_i * d_atm_ki(t) + atm_j * d_atm_kj(t) & - & + atm_k * d_atm_kk(t) + d_atm_extra_k(t)) * three_body_weight * rdecay(i,j) * rdecay(i,k) - & - & angular * radial(l) * rdecay(i,j) * d_ikdecay(t) * atm - - enddo - enddo - enddo - enddo - rep(i, twobody_size + 1:) = rep(i, twobody_size + 1:) + atom_rep - grad(i, twobody_size + 1:,:,:) = grad(i, twobody_size + 1:,:,:) + atom_grad -enddo + + dx = -(coordinates(i, k) - coordinates(j, k)) + + part(:) = ((log_Rs2(:) - mu)*(-dx*(rij**2*exp_s2 + eta2)/(rij*sqrt(exp_s2))**3) & + &*sqrt(exp_s2)/(sigma**2*rij) + (log_Rs2(:) - mu)**2*eta2*dx/ & + &(sigma**4*rij**4*exp_s2))*exp_ln/(Rs2(:)*sigma*sqrt(pi)*2) & + &- exp_ln*eta2*dx/(Rs2(:)*sigma**3*sqrt(pi)*rij**4*exp_s2*2.0d0) + + dscal = two_body_decay*dx/rij**(two_body_decay + 2.0d0) + ddecay = dx*0.5d0*pi*sin(pi*rij*invcut)*invcut*invrij + + part(:) = part(:)*scaling*rdecay(i, j) + radial_base(:)*dscal*rdecay(i, j) & + & + radial_base(:)*scaling*ddecay + + ! The gradients wrt coordinates + grad(i, (n - 1)*nbasis2 + 1:n*nbasis2, i, k) = grad(i, (n - 1)*nbasis2 + 1:n*nbasis2, i, k) + part + grad(i, (n - 1)*nbasis2 + 1:n*nbasis2, true_j, k) = grad(i, (n - 1)*nbasis2 + 1:n*nbasis2, true_j, k) - part + if (j <= natoms) then + grad(j, (m - 1)*nbasis2 + 1:m*nbasis2, i, k) = grad(j, (m - 1)*nbasis2 + 1:m*nbasis2, i, k) + part + add_grad(:, k, i, j) = part + end if + end do + end if + end do + end do !$OMP END PARALLEL DO -deallocate(rdecay) -deallocate(element_types) -deallocate(distance_matrix) -deallocate(inv_distance_matrix) -deallocate(sq_distance_matrix) -deallocate(inv_sq_distance_matrix) -deallocate(atom_rep) -deallocate(atom_grad) -deallocate(a) -deallocate(b) -deallocate(c) -deallocate(radial) -deallocate(angular_base) -deallocate(angular) -deallocate(d_angular) -deallocate(d_angular_d_i) -deallocate(d_angular_d_j) -deallocate(d_angular_d_k) -deallocate(d_radial) -deallocate(d_radial_d_i) -deallocate(d_radial_d_j) -deallocate(d_radial_d_k) -deallocate(d_ijdecay) -deallocate(d_ikdecay) +!$OMP PARALLEL DO PRIVATE(m) SCHEDULE(dynamic) + do j = 2, natoms + do i = 1, j - 1 + m = element_types(i) + rep(j, (m - 1)*nbasis2 + 1:m*nbasis2) = rep(j, (m - 1)*nbasis2 + 1:m*nbasis2) + add_rep(:, i, j) + do k = 1, 3 + grad(j, (m - 1)*nbasis2 + 1:m*nbasis2, j, k) = grad(j, (m - 1)*nbasis2 + 1:m*nbasis2, j, k) - add_grad(:, k, i, j) + end do + end do + end do +!$OMP END PARALLEL DO + + deallocate (add_rep, add_grad) + + deallocate (radial_base) + deallocate (radial) + deallocate (radial_part) + deallocate (part) + + ! Number of radial basis functions in the three body term + nbasis3 = size(Rs3) + ! Number of angular basis functions in the three body term + nabasis = size(Ts) + ! Size of two body terms + twobody_size = nelements*nbasis2 + + ! Inverse of the three body cutoff distance + invcut = 1.0d0/acut + ! Pre-calculate the three body decay + rdecay = decay(distance_matrix, invcut, natoms_tot) + + ! Allocate temporary + allocate (atom_rep(rep_size - twobody_size)) + allocate (atom_grad(rep_size - twobody_size, natoms, 3)) + allocate (a(3)) + allocate (b(3)) + allocate (c(3)) + allocate (radial(nbasis3)) + allocate (radial_base(nbasis3)) + allocate (angular_base(nabasis)) + allocate (angular(nabasis)) + allocate (d_angular(nabasis)) + allocate (d_angular_d_i(3)) + allocate (d_angular_d_j(3)) + allocate (d_angular_d_k(3)) + allocate (d_radial(nbasis3)) + allocate (d_radial_d_i(3)) + allocate (d_radial_d_j(3)) + allocate (d_radial_d_k(3)) + allocate (d_ijdecay(3)) + allocate (d_ikdecay(3)) + + allocate (d_atm_i(3)) + allocate (d_atm_j(3)) + allocate (d_atm_k(3)) + allocate (d_atm_ii(3)) + allocate (d_atm_ij(3)) + allocate (d_atm_ik(3)) + allocate (d_atm_ji(3)) + allocate (d_atm_jj(3)) + allocate (d_atm_jk(3)) + allocate (d_atm_ki(3)) + allocate (d_atm_kj(3)) + allocate (d_atm_kk(3)) + allocate (d_atm_i2(3)) + allocate (d_atm_j2(3)) + allocate (d_atm_k2(3)) + allocate (d_atm_i3(3)) + allocate (d_atm_j3(3)) + allocate (d_atm_k3(3)) + allocate (d_atm_extra_i(3)) + allocate (d_atm_extra_j(3)) + allocate (d_atm_extra_k(3)) + +!$OMP PARALLEL DO PRIVATE(atom_rep, atom_grad, rij, n, rij2, invrij, invrij2,& +!$OMP rik, m, rik2, invrik, invrjk, invrik2, a, b, c, angle, cos_i, cos_k,& +!$OMP cos_j, radial_base, radial, p, q, dot, angular, d_angular, d_angular_d_j,& +!$OMP d_angular_d_k, d_angular_d_i, d_radial, d_radial_d_j, d_radial_d_k,& +!$OMP d_radial_d_i, d_ijdecay, d_ikdecay, invr_atm, atm, atm_i, atm_j, atm_k, vi,& +!$OMP vj, vk, d_atm_ii, d_atm_ij, d_atm_ik, d_atm_ji, d_atm_jj, d_atm_jk, d_atm_ki,& +!$OMP d_atm_kj, d_atm_kk, d_atm_extra_i, d_atm_extra_j, d_atm_extra_k, s, z) SCHEDULE(dynamic) + do i = 1, natoms + atom_rep = 0.0d0 + atom_grad = 0.0d0 + do j = 1, natoms_tot - 1 + if (i .eq. j) cycle + ! distance between atom i and j + rij = distance_matrix(i, j) + if (rij > acut) cycle + ! index of the element of atom j + true_j = true_atom_id(j, natoms) + n = element_types(true_j) + ! squared distance between atom i and j + rij2 = sq_distance_matrix(i, j) + ! inverse distance between atom i and j + invrij = inv_distance_matrix(i, j) + ! inverse squared distance between atom i and j + invrij2 = inv_sq_distance_matrix(i, j) + do k = j + 1, natoms_tot + if (i .eq. k) cycle + ! distance between atom i and k + rik = distance_matrix(i, k) + if (rik > acut) cycle + ! index of the element of atom k + true_k = true_atom_id(k, natoms) + m = element_types(true_k) + ! squared distance between atom i and k + rik2 = sq_distance_matrix(i, k) + ! inverse distance between atom i and k + invrik = inv_distance_matrix(i, k) + ! inverse distance between atom j and k + invrjk = inv_distance_matrix(j, k) + ! inverse squared distance between atom i and k + invrik2 = inv_sq_distance_matrix(i, k) + ! coordinates of atoms j, i, k + a = coordinates(j, :) + b = coordinates(i, :) + c = coordinates(k, :) + ! angle between atom i, j and k, centered on i + angle = calc_angle(a, b, c) + cos_i = calc_cos_angle(a, b, c) + cos_k = calc_cos_angle(a, c, b) + cos_j = calc_cos_angle(b, a, c) + + ! part of the radial part of the 3body terms + radial_base(:) = exp(-eta3*(0.5d0*(rij + rik) - Rs3(:))**2) + radial(:) = radial_base(:) ! * scaling + + p = min(n, m) - 1 + ! the highest index of the elements of j,k + q = max(n, m) - 1 + ! Dot product between the vectors connecting atom i,j and i,k + dot = dot_product(a - b, c - b) + + angular(1) = exp(-(zeta**2)*0.5d0)*2*cos(angle) + angular(2) = exp(-(zeta**2)*0.5d0)*2*sin(angle) + + d_angular(1) = exp(-(zeta**2)*0.5d0)*2*sin(angle)/sqrt(max(1d-10, rij2*rik2 - dot**2)) + d_angular(2) = -exp(-(zeta**2)*0.5d0)*2*cos(angle)/sqrt(max(1d-10, rij2*rik2 - dot**2)) + + ! Part of the derivative of the angular basis functions wrt atom j (dim(3)) + d_angular_d_j = c - b + dot*((b - a)*invrij2) + ! Part of the derivative of the angular basis functions wrt atom k (dim(3)) + d_angular_d_k = a - b + dot*((b - c)*invrik2) + ! Part of the derivative of the angular basis functions wrt atom i (dim(3)) + d_angular_d_i = -(d_angular_d_j + d_angular_d_k) + + ! Part of the derivative of the radial basis functions wrt coordinates (dim(nbasis3)) + ! including decay + d_radial = radial*eta3*(0.5d0*(rij + rik) - Rs3) ! * rdecay(i,j) * rdecay(i,k) + ! Part of the derivative of the radial basis functions wrt atom j (dim(3)) + d_radial_d_j = (b - a)*invrij + ! Part of the derivative of the radial basis functions wrt atom k (dim(3)) + d_radial_d_k = (b - c)*invrik + ! Part of the derivative of the radial basis functions wrt atom i (dim(3)) + d_radial_d_i = -(d_radial_d_j + d_radial_d_k) + + ! Part of the derivative of the i,j decay functions wrt coordinates (dim(3)) + d_ijdecay = -pi*(b - a)*sin(pi*rij*invcut)*0.5d0*invrij*invcut + ! Part of the derivative of the i,k decay functions wrt coordinates (dim(3)) + d_ikdecay = -pi*(b - c)*sin(pi*rik*invcut)*0.5d0*invrik*invcut + + invr_atm = (invrij*invrjk*invrik)**three_body_decay + + ! Axilrod-Teller-Muto term + atm = (1.0d0 + 3.0d0*cos_i*cos_j*cos_k)*invr_atm*three_body_weight + + atm_i = (3.0d0*cos_j*cos_k)*invr_atm*invrij*invrik + atm_j = (3.0d0*cos_k*cos_i)*invr_atm*invrij*invrjk + atm_k = (3.0d0*cos_i*cos_j)*invr_atm*invrjk*invrik + + vi = dot_product(a - b, c - b) + vj = dot_product(c - a, b - a) + vk = dot_product(b - c, a - c) + + d_atm_ii(:) = 2*b - a - c - vi*((b - a)*invrij**2 + (b - c)*invrik**2) + d_atm_ij(:) = c - a - vj*(b - a)*invrij**2 + d_atm_ik(:) = a - c - vk*(b - c)*invrik**2 + + d_atm_ji(:) = c - b - vi*(a - b)*invrij**2 + d_atm_jj(:) = 2*a - b - c - vj*((a - b)*invrij**2 + (a - c)*invrjk**2) + d_atm_jk(:) = b - c - vk*(a - c)*invrjk**2 + + d_atm_ki(:) = a - b - vi*(c - b)*invrik**2 + d_atm_kj(:) = b - a - vj*(c - a)*invrjk**2 + d_atm_kk(:) = 2*c - a - b - vk*((c - a)*invrjk**2 + (c - b)*invrik**2) + + d_atm_extra_i(:) = ((a - b)*invrij**2 + (c - b)*invrik**2)*atm*three_body_decay/three_body_weight + d_atm_extra_j(:) = ((b - a)*invrij**2 + (c - a)*invrjk**2)*atm*three_body_decay/three_body_weight + d_atm_extra_k(:) = ((a - c)*invrjk**2 + (b - c)*invrik**2)*atm*three_body_decay/three_body_weight + + ! Get index of where the contributions of atoms i,j,k should be added + s = nbasis3*nabasis*(-(p*(p + 1))/2 + q + nelements*p) + 1 + + do l = 1, nbasis3 + + ! Get index of where the contributions of atoms i,j,k should be added + z = s + (l - 1)*nabasis + + ! Add the contributions for atoms i,j,k + atom_rep(z:z + nabasis - 1) = atom_rep(z:z + nabasis - 1) & + & + angular*radial(l)*atm*rdecay(i, j)*rdecay(i, k) + + do t = 1, 3 + + ! Add up all gradient contributions wrt atom i + atom_grad(z:z + nabasis - 1, i, t) = atom_grad(z:z + nabasis - 1, i, t) + & + & d_angular*d_angular_d_i(t)*radial(l)*atm*rdecay(i, j)*rdecay(i, k) + & + & angular*d_radial(l)*d_radial_d_i(t)*atm*rdecay(i, j)*rdecay(i, k) + & + & angular*radial(l)*(atm_i*d_atm_ii(t) + atm_j*d_atm_ij(t) & + & + atm_k*d_atm_ik(t) + d_atm_extra_i(t))*three_body_weight*rdecay(i, j)*rdecay(i, k) + & + & angular*radial(l)*(d_ijdecay(t)*rdecay(i, k) + rdecay(i, j)*d_ikdecay(t))*atm + ! Add up all gradient contributions wrt atom j + atom_grad(z:z + nabasis - 1, true_j, t) = atom_grad(z:z + nabasis - 1, true_j, t) + & + & d_angular*d_angular_d_j(t)*radial(l)*atm*rdecay(i, j)*rdecay(i, k) + & + & angular*d_radial(l)*d_radial_d_j(t)*atm*rdecay(i, j)*rdecay(i, k) + & + & angular*radial(l)*(atm_i*d_atm_ji(t) + atm_j*d_atm_jj(t) & + & + atm_k*d_atm_jk(t) + d_atm_extra_j(t))*three_body_weight*rdecay(i, j)*rdecay(i, k) - & + & angular*radial(l)*d_ijdecay(t)*rdecay(i, k)*atm + + ! Add up all gradient contributions wrt atom k + atom_grad(z:z + nabasis - 1, true_k, t) = atom_grad(z:z + nabasis - 1, true_k, t) + & + & d_angular*d_angular_d_k(t)*radial(l)*atm*rdecay(i, j)*rdecay(i, k) + & + & angular*d_radial(l)*d_radial_d_k(t)*atm*rdecay(i, j)*rdecay(i, k) + & + & angular*radial(l)*(atm_i*d_atm_ki(t) + atm_j*d_atm_kj(t) & + & + atm_k*d_atm_kk(t) + d_atm_extra_k(t))*three_body_weight*rdecay(i, j)*rdecay(i, k) - & + & angular*radial(l)*rdecay(i, j)*d_ikdecay(t)*atm + + end do + end do + end do + end do + rep(i, twobody_size + 1:) = rep(i, twobody_size + 1:) + atom_rep + grad(i, twobody_size + 1:, :, :) = grad(i, twobody_size + 1:, :, :) + atom_grad + end do +!$OMP END PARALLEL DO + deallocate (rdecay) + deallocate (element_types) + deallocate (distance_matrix) + deallocate (inv_distance_matrix) + deallocate (sq_distance_matrix) + deallocate (inv_sq_distance_matrix) + deallocate (atom_rep) + deallocate (atom_grad) + deallocate (a) + deallocate (b) + deallocate (c) + deallocate (radial) + deallocate (angular_base) + deallocate (angular) + deallocate (d_angular) + deallocate (d_angular_d_i) + deallocate (d_angular_d_j) + deallocate (d_angular_d_k) + deallocate (d_radial) + deallocate (d_radial_d_i) + deallocate (d_radial_d_j) + deallocate (d_radial_d_k) + deallocate (d_ijdecay) + deallocate (d_ikdecay) end subroutine fgenerate_fchl_acsf_and_gradients