diff --git a/qml/kernels/fgradient_kernels.f90 b/qml/kernels/fgradient_kernels.f90 index a65f8a57b..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 @@ -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..7d60efd0a 100644 --- a/qml/representations/facsf.f90 +++ b/qml/representations/facsf.f90 @@ -1,292 +1,278 @@ -! MIT License -! -! Copyright (c) 2018 Lars Andersen Bratholm -! -! 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 - 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 -end module acsf_utils + if (atom_id <= natoms) then + tatom_id = atom_id + else + tatom_id = modulo(atom_id - 1, natoms) + 1 + end if + + 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 @@ -294,1012 +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, rep_size, & - & two_body_decay, three_body_decay, three_body_weight, rep) - - use acsf_utils, only: decay, calc_angle, calc_cos_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(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 :: 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 :: 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)) - - ! 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 - continue - 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)) - - radial = 0.0d0 - ! !$OMP PARALLEL DO PRIVATE(n,m,rij,radial) REDUCTION(+:rep) - do i = 1, natoms - ! index of the element of atom i - m = element_types(i) - do j = i + 1, natoms - ! 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 - - 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 - endif - enddo - enddo - ! !$OMP END PARALLEL DO - - deallocate(radial) - - ! 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(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) - 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 - 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 = 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 - rep3(i, z:z + nabasis - 1) = rep3(i, z:z + nabasis - 1) + angular * radial(l) * ksi3 - enddo - enddo - enddo - enddo - ! !$OMP END PARALLEL DO - - rep = rep + rep3 - - deallocate(element_types) - deallocate(rdecay) - deallocate(distance_matrix) - deallocate(rep3) - deallocate(a) - deallocate(b) - deallocate(c) - deallocate(radial) - deallocate(angular) + & 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)) + 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 + 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, 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 - - 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 - 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, 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 - 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 - do i = 1, natoms - do j = 1, nelements - if (nuclear_charges(i) .eq. elements(j)) then - element_types(i) = j - continue - 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 = i+1, natoms - 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(exp_ln(nbasis2)) - allocate(log_Rs2(nbasis2)) - - grad =0.0d0 - rep = 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) - 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) - - 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 - - - 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, 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 - - enddo - endif - enddo - enddo - ! !$OMP END PARALLEL DO - - 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) - - ! 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)) - - ! ! 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 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, 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) + & - & 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, 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) + & - & 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 - ! !$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, 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 + 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) + + 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 + 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 + + 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 + +!$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 diff --git a/qml/representations/representations.py b/qml/representations/representations.py index 9cf89c360..35610d68b 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 @@ -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(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): + 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) 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 ...") diff --git a/test/test_fchl_acsf_pbc.py b/test/test_fchl_acsf_pbc.py new file mode 100644 index 000000000..64c52b983 --- /dev/null +++ b/test/test_fchl_acsf_pbc.py @@ -0,0 +1,136 @@ +# 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 +REP_PARAMS["acut"]=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 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) + all_charges=deepcopy(nuclear_charges) + 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): + 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= pbc_corrected_drep(drep, 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, 16) + + 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() +