Skip to content

Commit

Permalink
Add analytical stress tensor
Browse files Browse the repository at this point in the history
  • Loading branch information
@cschran @dev-zero
cschran authored and dev-zero committed Dec 8, 2021
1 parent e5dc313 commit b5b8227
Show file tree
Hide file tree
Showing 6 changed files with 383 additions and 9 deletions.
63 changes: 58 additions & 5 deletions src/nnp_acsf.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -51,15 +51,16 @@ MODULE nnp_acsf
!> \param nnp ...
!> \param i ...
!> \param dsymdxyz ...
!> \param stress ...
!> \date 2020-10-10
!> \author Christoph Schran (christoph.schran@rub.de)
! **************************************************************************************************
SUBROUTINE nnp_calc_acsf(nnp, i, dsymdxyz)
SUBROUTINE nnp_calc_acsf(nnp, i, dsymdxyz, stress)

TYPE(nnp_type), INTENT(INOUT), POINTER :: nnp
INTEGER, INTENT(IN) :: i
REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT), &
OPTIONAL :: dsymdxyz
OPTIONAL :: dsymdxyz, stress

CHARACTER(len=*), PARAMETER :: routineN = 'nnp_calc_acsf', routineP = moduleN//':'//routineN

Expand Down Expand Up @@ -105,6 +106,11 @@ SUBROUTINE nnp_calc_acsf(nnp, i, dsymdxyz)
dsymdxyz(l, m, i) = dsymdxyz(l, m, i) + forcetmp(l, sf)
dsymdxyz(l, m, jj) = dsymdxyz(l, m, jj) - forcetmp(l, sf)
END DO
IF (PRESENT(stress)) THEN
DO l = 1, 3
stress(:, l, m) = stress(:, l, m) + rvect1(:)*forcetmp(l, sf)
END DO
END IF
nnp%rad(ind)%y(m) = nnp%rad(ind)%y(m) + symtmp(sf)
END DO
END DO
Expand Down Expand Up @@ -144,6 +150,12 @@ SUBROUTINE nnp_calc_acsf(nnp, i, dsymdxyz)
dsymdxyz(l, m, kk) = dsymdxyz(l, m, kk) &
+ force3tmp(l, 3, sf)
END DO
IF (PRESENT(stress)) THEN
DO l = 1, 3
stress(:, l, m) = stress(:, l, m) - rvect1(:)*force3tmp(l, 2, sf)
stress(:, l, m) = stress(:, l, m) - rvect2(:)*force3tmp(l, 3, sf)
END DO
END IF
nnp%ang(ind)%y(m - off) = nnp%ang(ind)%y(m - off) + symtmp(sf)
END DO
END IF
Expand Down Expand Up @@ -177,6 +189,12 @@ SUBROUTINE nnp_calc_acsf(nnp, i, dsymdxyz)
dsymdxyz(l, m, kk) = dsymdxyz(l, m, kk) &
+ force3tmp(l, 3, sf)
END DO
IF (PRESENT(stress)) THEN
DO l = 1, 3
stress(:, l, m) = stress(:, l, m) - rvect1(:)*force3tmp(l, 2, sf)
stress(:, l, m) = stress(:, l, m) - rvect2(:)*force3tmp(l, 3, sf)
END DO
END IF
nnp%ang(ind)%y(m - off) = nnp%ang(ind)%y(m - off) + symtmp(sf)
END DO
END IF
Expand Down Expand Up @@ -259,7 +277,11 @@ SUBROUTINE nnp_calc_acsf(nnp, i, dsymdxyz)
CALL nnp_check_extrapolation(nnp, ind)

IF (PRESENT(dsymdxyz)) THEN
CALL nnp_scale_acsf(nnp, ind, dsymdxyz)
IF (PRESENT(stress)) THEN
CALL nnp_scale_acsf(nnp, ind, dsymdxyz, stress)
ELSE
CALL nnp_scale_acsf(nnp, ind, dsymdxyz)
END IF
ELSE
CALL nnp_scale_acsf(nnp, ind)
END IF
Expand Down Expand Up @@ -316,15 +338,16 @@ END SUBROUTINE nnp_check_extrapolation
!> \param nnp ...
!> \param ind ...
!> \param dsymdxyz ...
!> \param stress ...
!> \date 2020-10-10
!> \author Christoph Schran (christoph.schran@rub.de)
! **************************************************************************************************
SUBROUTINE nnp_scale_acsf(nnp, ind, dsymdxyz)
SUBROUTINE nnp_scale_acsf(nnp, ind, dsymdxyz, stress)

TYPE(nnp_type), INTENT(INOUT) :: nnp
INTEGER, INTENT(IN) :: ind
REAL(KIND=dp), DIMENSION(:, :, :), INTENT(OUT), &
OPTIONAL :: dsymdxyz
OPTIONAL :: dsymdxyz, stress

INTEGER :: j, k, off

Expand Down Expand Up @@ -430,6 +453,36 @@ SUBROUTINE nnp_scale_acsf(nnp, ind, dsymdxyz)
END IF
END IF

IF (PRESENT(stress)) THEN
IF (nnp%scale_acsf) THEN
DO j = 1, nnp%n_rad(ind)
stress(:, :, j) = stress(:, :, j)/ &
(nnp%rad(ind)%loc_max(j) - &
nnp%rad(ind)%loc_min(j))* &
(nnp%scmax - nnp%scmin)
END DO
off = nnp%n_rad(ind)
DO j = 1, nnp%n_ang(ind)
stress(:, :, j + off) = stress(:, :, j + off)/ &
(nnp%ang(ind)%loc_max(j) - &
nnp%ang(ind)%loc_min(j))* &
(nnp%scmax - nnp%scmin)
END DO
ELSE IF (nnp%scale_sigma_acsf) THEN
DO j = 1, nnp%n_rad(ind)
stress(:, :, j) = stress(:, :, j)/ &
nnp%rad(ind)%sigma(j)* &
(nnp%scmax - nnp%scmin)
END DO
off = nnp%n_rad(ind)
DO j = 1, nnp%n_ang(ind)
stress(:, :, j + off) = stress(:, :, j + off)/ &
nnp%ang(ind)%sigma(j)* &
(nnp%scmax - nnp%scmin)
END DO
END IF
END IF

END SUBROUTINE nnp_scale_acsf

! **************************************************************************************************
Expand Down
1 change: 1 addition & 0 deletions src/nnp_environment.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -287,6 +287,7 @@ SUBROUTINE nnp_init_model(nnp_env)
ALLOCATE (nnp_env%committee_energy(nnp_env%n_committee))
ALLOCATE (nnp_env%myforce(3, nnp_env%num_atoms, nnp_env%n_committee))
ALLOCATE (nnp_env%committee_forces(3, nnp_env%num_atoms, nnp_env%n_committee))
ALLOCATE (nnp_env%committee_stress(3, 3, nnp_env%n_committee))

CALL section_vals_val_get(nnp_env%nnp_input, "NNP_INPUT_FILE_NAME", c_val=file_name)
CALL parser_create(parser, file_name, para_env=logger%para_env)
Expand Down
3 changes: 2 additions & 1 deletion src/nnp_environment_types.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -119,7 +119,7 @@ MODULE nnp_environment_types
INTEGER, ALLOCATABLE, DIMENSION(:) :: ele_ind, nuc_atoms, sort, sort_inv
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: coord
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: myforce
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: committee_forces
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: committee_forces, committee_stress
CHARACTER(len=default_string_length), &
ALLOCATABLE, DIMENSION(:) :: atoms
REAL(KIND=dp), DIMENSION(:, :), ALLOCATABLE :: nnp_forces
Expand Down Expand Up @@ -413,6 +413,7 @@ SUBROUTINE nnp_env_release(nnp_env)
IF (ALLOCATED(nnp_env%coord)) DEALLOCATE (nnp_env%coord)
IF (ALLOCATED(nnp_env%myforce)) DEALLOCATE (nnp_env%myforce)
IF (ALLOCATED(nnp_env%committee_forces)) DEALLOCATE (nnp_env%committee_forces)
IF (ALLOCATED(nnp_env%committee_stress)) DEALLOCATE (nnp_env%committee_stress)
IF (ALLOCATED(nnp_env%atoms)) DEALLOCATE (nnp_env%atoms)
IF (ALLOCATED(nnp_env%nnp_forces)) DEALLOCATE (nnp_env%nnp_forces)

Expand Down
35 changes: 32 additions & 3 deletions src/nnp_force.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,7 +19,8 @@ MODULE nnp_force
cp_print_key_finished_output,&
cp_print_key_should_output,&
cp_print_key_unit_nr
USE cp_subsys_types, ONLY: cp_subsys_type
USE cp_subsys_types, ONLY: cp_subsys_get,&
cp_subsys_type
USE cp_units, ONLY: cp_unit_from_cp2k
USE distribution_1d_types, ONLY: distribution_1d_type
USE input_section_types, ONLY: section_vals_get_subs_vals,&
Expand All @@ -38,6 +39,7 @@ MODULE nnp_force
USE particle_types, ONLY: particle_type
USE periodic_table, ONLY: get_ptable_info
USE physcon, ONLY: angstrom
USE virial_types, ONLY: virial_type
#include "./base/base_uses.f90"

IMPLICIT NONE
Expand Down Expand Up @@ -68,14 +70,16 @@ SUBROUTINE nnp_calc_energy_force(nnp, calc_forces)
INTEGER :: handle, i, i_com, ig, ind, istart, j, k, &
m, mecalc
INTEGER, ALLOCATABLE, DIMENSION(:) :: allcalc
LOGICAL :: calc_stress
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: denergydsym
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: dsymdxyz
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: dsymdxyz, stress
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
TYPE(cp_logger_type), POINTER :: logger
TYPE(cp_subsys_type), POINTER :: subsys
TYPE(distribution_1d_type), POINTER :: local_particles
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
TYPE(section_vals_type), POINTER :: print_section
TYPE(virial_type), POINTER :: virial

CALL timeset(routineN, handle)

Expand All @@ -88,10 +92,17 @@ SUBROUTINE nnp_calc_energy_force(nnp, calc_forces)
subsys=subsys, local_particles=local_particles, &
atomic_kind_set=atomic_kind_set)

CALL cp_subsys_get(subsys, &
virial=virial)

calc_stress = virial%pv_availability .AND. (.NOT. virial%pv_numer)
IF (calc_stress .AND. .NOT. calc_forces) CPABORT('Stress cannot be calculated without forces')

! Initialize energy and gradient:
nnp%atomic_energy(:, :) = 0.0_dp
IF (calc_forces) nnp%myforce(:, :, :) = 0.0_dp
IF (calc_forces) nnp%committee_forces(:, :, :) = 0.0_dp
IF (calc_stress) nnp%committee_stress(:, :, :) = 0.0_dp

!fill coord array
ig = 1
Expand Down Expand Up @@ -138,7 +149,13 @@ SUBROUTINE nnp_calc_energy_force(nnp, calc_forces)
nnp%arc(ind)%layer(1)%node_grad(:) = 0.0_dp
ALLOCATE (dsymdxyz(3, nnp%arc(ind)%n_nodes(1), nnp%num_atoms))
dsymdxyz(:, :, :) = 0.0_dp
CALL nnp_calc_acsf(nnp, i, dsymdxyz)
IF (calc_stress) THEN
ALLOCATE (stress(3, 3, nnp%arc(ind)%n_nodes(1)))
stress(:, :, :) = 0.0_dp
CALL nnp_calc_acsf(nnp, i, dsymdxyz, stress)
ELSE
CALL nnp_calc_acsf(nnp, i, dsymdxyz)
END IF
ELSE
CALL nnp_calc_acsf(nnp, i)
END IF
Expand All @@ -160,6 +177,10 @@ SUBROUTINE nnp_calc_energy_force(nnp, calc_forces)
nnp%myforce(m, k, i_com) = nnp%myforce(m, k, i_com) - denergydsym(j)*dsymdxyz(m, j, k)
END DO
END DO
IF (calc_stress) THEN
nnp%committee_stress(:, :, i_com) = nnp%committee_stress(:, :, i_com) - &
denergydsym(j)*stress(:, :, j)
END IF
END DO
DEALLOCATE (denergydsym)
END IF
Expand All @@ -168,6 +189,9 @@ SUBROUTINE nnp_calc_energy_force(nnp, calc_forces)
!deallocate memory
IF (calc_forces) THEN
DEALLOCATE (dsymdxyz)
IF (calc_stress) THEN
DEALLOCATE (stress)
END IF
END IF

END DO ! loop over num_atoms
Expand All @@ -191,6 +215,11 @@ SUBROUTINE nnp_calc_energy_force(nnp, calc_forces)
END DO
END IF

IF (calc_stress) THEN
CALL mp_sum(nnp%committee_stress, logger%para_env%group)
virial%pv_virial = SUM(nnp%committee_stress, DIM=3)/REAL(nnp%n_committee, dp)
END IF

! Bias the standard deviation of committee disagreement
IF (nnp%bias) THEN
CALL nnp_bias_sigma(nnp, calc_forces)
Expand Down

0 comments on commit b5b8227

Please sign in to comment.