From e22cb24e4437376446a0940e2384ad4923423b54 Mon Sep 17 00:00:00 2001 From: Albert Bartok-Partay Date: Sun, 5 Jul 2026 15:15:48 +0100 Subject: [PATCH 1/2] Dipole and local dipole prediction added. --- src/Potentials/IPModel_GAP.F90 | 65 +++++++++++++++++++---- src/Potentials/Potential_Sum_routines.F90 | 36 +++++++++++-- src/Potentials/Potential_simple.F90 | 23 ++++++-- 3 files changed, 107 insertions(+), 17 deletions(-) diff --git a/src/Potentials/IPModel_GAP.F90 b/src/Potentials/IPModel_GAP.F90 index 1352d5c6dd..415cc82e5b 100644 --- a/src/Potentials/IPModel_GAP.F90 +++ b/src/Potentials/IPModel_GAP.F90 @@ -249,22 +249,25 @@ subroutine IPModel_GAP_Calc(this, at, e, local_e, f, virial, local_virial, args_ real(dp), dimension(:,:), allocatable :: f_in - real(dp), dimension(3) :: pos, f_gp + real(dp), dimension(3) :: pos, f_gp, dipole, dipole_i real(dp), dimension(3,3) :: virial_i type(Dictionary) :: params logical, dimension(:), pointer :: atom_mask_pointer logical, dimension(:), allocatable :: mpi_local_mask logical :: has_atom_mask_name - character(STRING_LENGTH) :: atom_mask_name, calc_local_gap_variance, calc_energy_per_coordinate + character(STRING_LENGTH) :: atom_mask_name, calc_local_gap_variance, calc_energy_per_coordinate, & + calc_dipole, calc_local_dipole real(dp) :: r_scale, E_scale real(dp) :: gap_variance_i_cutoff real(dp), dimension(:), allocatable :: gap_variance, local_gap_variance_in real(dp), dimension(:), pointer :: local_gap_variance_pointer real(dp), dimension(:,:), allocatable :: gap_variance_gradient_in - real(dp), dimension(:,:), pointer :: gap_variance_gradient_pointer + real(dp), dimension(:,:), pointer :: gap_variance_gradient_pointer, at_local_dipole_ptr real(dp) :: gap_variance_regularisation - logical :: do_rescale_r, do_rescale_E, do_gap_variance, print_gap_variance, do_local_gap_variance, do_energy_per_coordinate + logical :: do_rescale_r, do_rescale_E, do_gap_variance, print_gap_variance, & + do_local_gap_variance, do_energy_per_coordinate, do_grad_descriptor + logical :: do_calc_dipole, do_calc_local_dipole integer :: only_descriptor logical :: do_select_descriptor logical :: mpi_parallel_descriptor @@ -319,6 +322,9 @@ subroutine IPModel_GAP_Calc(this, at, e, local_e, f, virial, local_virial, args_ atom_mask_name = "" only_descriptor = 0 + calc_dipole = "" + calc_local_dipole = "" + call initialise(params) call param_register(params, 'atom_mask_name', 'NONE',atom_mask_name,has_value_target=has_atom_mask_name, & @@ -336,6 +342,10 @@ subroutine IPModel_GAP_Calc(this, at, e, local_e, f, virial, local_virial, args_ call param_register(params, 'mpi_parallel_descriptor', 'F', mpi_parallel_descriptor, help_string="Do MPI parallelism over descriptor instances rather than atoms") + call param_register(params, "dipole", "", calc_dipole, help_string="If present, calculate dipole and put it in field with this string as name") + + call param_register(params, "local_dipole", "", calc_local_dipole, help_string="If present, calculate local dipole and put it in field with this string as name") + if(present(args_str)) then if (.not. param_read_line(params,args_str,ignore_unknown=.true.,task='IPModel_GAP_Calc args_str')) & call system_abort("IPModel_GAP_Calc failed to parse args_str='"//trim(args_str)//"'") @@ -390,6 +400,14 @@ subroutine IPModel_GAP_Calc(this, at, e, local_e, f, virial, local_virial, args_ endif endif + do_calc_dipole = len_trim(calc_dipole) > 0 + do_calc_local_dipole = len_trim(calc_local_dipole) > 0 + + if( do_calc_local_dipole ) then + call add_property(at, trim(calc_local_dipole), 0.0_dp, n_cols=3, ptr2=at_local_dipole_ptr) + at_local_dipole_ptr = 0.0_dp + endif + dipole = 0.0_dp if(print_gap_variance) then call print('GAP_VARIANCE potential '//trim(this%label)//' calculating for '//this%my_gp%n_coordinate//' descriptors') @@ -413,23 +431,29 @@ subroutine IPModel_GAP_Calc(this, at, e, local_e, f, virial, local_virial, args_ call gpCoordinates_initialise_variance_estimate(this%my_gp%coordinate(i_coordinate), gap_variance_regularisation) endif - if(present(f) .or. present(virial) .or. present(local_virial)) then + do_grad_descriptor = present(f) .or. present(virial) .or. present(local_virial) .or. & + do_calc_dipole .or. do_calc_local_dipole + + if( do_grad_descriptor ) then if (allocated(gradPredict)) deallocate(gradPredict) allocate(gradPredict(d)) if(allocated(grad_variance_estimate)) deallocate(grad_variance_estimate) allocate(grad_variance_estimate(d)) end if + call calc(this%my_descriptor(i_coordinate),at,my_descriptor_data, & - do_descriptor=.true.,do_grad_descriptor=present(f) .or. present(virial) .or. present(local_virial), args_str=trim(string(my_args_str)), error=error) + do_descriptor=.true.,do_grad_descriptor=do_grad_descriptor, & + args_str=trim(string(my_args_str)), error=error) + PASS_ERROR(error) allocate(gap_variance(size(my_descriptor_data%x))) call system_timer('IPModel_GAP_Calc_gp_predict') !$omp parallel default(none) private(i,gradPredict, grad_variance_estimate, e_i,n,m,j,pos,f_gp,e_i_cutoff,virial_i,i_pos0,gap_variance_i_cutoff) & -!$omp shared(this,at,i_coordinate,my_descriptor_data,e,virial,local_virial,local_e,do_gap_variance,do_local_gap_variance,gap_variance,f,do_energy_per_coordinate,mpi,mpi_parallel_descriptor) & -!$omp reduction(+:local_e_in,f_in,virial_in,local_gap_variance_in, gap_variance_gradient_in, energy_per_coordinate) +!$omp shared(this,at,i_coordinate,my_descriptor_data,e,virial,local_virial,local_e,do_gap_variance,do_local_gap_variance,gap_variance,f,do_energy_per_coordinate,mpi,mpi_parallel_descriptor,at_local_energy_ptr) & +!$omp reduction(+:local_e_in,f_in,virial_in,local_gap_variance_in, gap_variance_gradient_in, energy_per_coordinate,dipole) !$omp do schedule(dynamic) loop_over_descriptor_instances: do i = 1, size(my_descriptor_data%x) @@ -442,7 +466,7 @@ subroutine IPModel_GAP_Calc(this, at, e, local_e, f, virial, local_virial, args_ !call system_timer('IPModel_GAP_Calc_gp_predict') - if(present(f) .or. present(virial) .or. present(local_virial)) then + if( do_grad_descriptor ) then call reallocate(gradPredict,size(my_descriptor_data%x(i)%data(:)),zero=.true.) e_i = gp_predict(this%my_gp%coordinate(i_coordinate) , xStar=my_descriptor_data%x(i)%data(:), gradPredict = gradPredict, variance_estimate=gap_variance(i), do_variance_estimate=do_gap_variance, grad_variance_estimate=grad_variance_estimate) else @@ -469,6 +493,28 @@ subroutine IPModel_GAP_Calc(this, at, e, local_e, f, virial, local_virial, args_ enddo endif + if( do_calc_dipole .or. do_calc_local_dipole ) then + if( size( my_descriptor_data%x(i)%ci ) /= 1 ) then + RAISE_ERROR("IPModel_GAP_Calc: descriptor must have a well-defined central atom, try e.g. SOAP",error) + endif + + j = my_descriptor_data%x(i)%ci(1) + dipole_i = 0.0_dp + + do n = lbound(my_descriptor_data%x(i)%ii,1), ubound(my_descriptor_data%x(i)%ii,1) + if( .not. my_descriptor_data%x(i)%has_grad_data(n) .or. & + my_descriptor_data%x(i)%ci(1) /= my_descriptor_data%x(i)%ii(n) ) cycle + + dipole_i = matmul( gradPredict,my_descriptor_data%x(i)%grad_data(:,:,n)) * my_descriptor_data%x(i)%covariance_cutoff + & + e_i * my_descriptor_data%x(i)%grad_covariance_cutoff(:,n) + + enddo + if( do_calc_dipole ) dipole = dipole + dipole_i + if( do_calc_local_dipole ) at_local_dipole_ptr(:,my_descriptor_data%x(i)%ci(1)) = & + at_local_dipole_ptr(:,my_descriptor_data%x(i)%ci(1)) + dipole_i + + endif + if(present(f) .or. present(virial) .or. present(local_virial)) then i_pos0 = lbound(my_descriptor_data%x(i)%ii,1) @@ -524,6 +570,7 @@ subroutine IPModel_GAP_Calc(this, at, e, local_e, f, virial, local_virial, args_ if(present(e)) e = sum(local_e_in) if(present(local_e)) local_e = local_e_in if(present(virial)) virial = sum(virial_in,dim=3) + if(do_calc_dipole) call set_value(at%params, trim(calc_dipole), dipole) if(present(local_virial)) then do i = 1, at%N diff --git a/src/Potentials/Potential_Sum_routines.F90 b/src/Potentials/Potential_Sum_routines.F90 index bd2f3646ca..4db57cfced 100644 --- a/src/Potentials/Potential_Sum_routines.F90 +++ b/src/Potentials/Potential_Sum_routines.F90 @@ -63,15 +63,16 @@ recursive subroutine Potential_Sum_Calc(this, at, args_str, error) character(*), intent(in), optional :: args_str integer, intent(out), optional :: error - real(dp) :: energy, virial(3,3) - real(dp), pointer :: at_force_ptr(:,:), at_local_energy_ptr(:), at_local_virial_ptr(:,:) + real(dp) :: energy, virial(3,3), dipole(3) + real(dp), pointer :: at_force_ptr(:,:), at_local_energy_ptr(:), at_local_virial_ptr(:,:), at_local_dipole_ptr(:,:) real(dp) :: my_e_1, my_e_2 real(dp), allocatable :: my_local_e_1(:) - real(dp), allocatable :: my_f_1(:,:), my_local_virial_1(:,:) - real(dp) :: my_virial_1(3,3) + real(dp), allocatable :: my_f_1(:,:), my_local_virial_1(:,:), my_local_dipole_1(:,:) + real(dp) :: my_virial_1(3,3), my_dipole_1(3) type(Dictionary) :: params - character(STRING_LENGTH) :: calc_energy, calc_force, calc_local_energy, calc_virial, calc_local_virial, calc_args_pot1, calc_args_pot2, my_args_str + character(STRING_LENGTH) :: calc_energy, calc_force, calc_local_energy, calc_virial, calc_local_virial, & + calc_dipole, calc_local_dipole, calc_args_pot1, calc_args_pot2, my_args_str logical :: store_contributions INIT_ERROR(error) @@ -80,8 +81,10 @@ recursive subroutine Potential_Sum_Calc(this, at, args_str, error) call param_register(params,"energy", "", calc_energy, help_string="No help yet. This source file was $LastChangedBy$") call param_register(params,"force", "", calc_force, help_string="No help yet. This source file was $LastChangedBy$") call param_register(params,"virial", "", calc_virial, help_string="No help yet. This source file was $LastChangedBy$") + call param_register(params,"dipole", "", calc_dipole, help_string="Whether to force calculation of the dipole.") call param_register(params,"local_energy", "", calc_local_energy, help_string="No help yet. This source file was $LastChangedBy$") call param_register(params,"local_virial", "", calc_local_virial, help_string="No help yet. This source file was $LastChangedBy$") + call param_register(params,"local_dipole", "", calc_local_dipole, help_string="Whether to force calculation of the local dipole.") call param_register(params,"calc_args_pot1", "", calc_args_pot1, help_string="additional args_str to pass along to pot1") call param_register(params,"calc_args_pot2", "", calc_args_pot2, help_string="additional args_str to pass along to pot2") call param_register(params,"store_contributions", "F", store_contributions, help_string="if true, store contributions to sum with _pot1 and _pot2 suffixes") @@ -103,6 +106,10 @@ recursive subroutine Potential_Sum_Calc(this, at, args_str, error) call get_param_value(at, trim(calc_virial), my_virial_1) if (store_contributions) call set_param_value(at, trim(calc_virial)//"_pot1", my_virial_1) endif + if (len_trim(calc_dipole) > 0) then + call get_param_value(at, trim(calc_dipole), my_dipole_1) + if (store_contributions) call set_param_value(at, trim(calc_dipole)//"_pot1", my_dipole_1) + endif if (len_trim(calc_local_energy) > 0) then call assign_property_pointer(at, trim(calc_local_energy), at_local_energy_ptr, error=error) PASS_ERROR(error) @@ -124,6 +131,13 @@ recursive subroutine Potential_Sum_Calc(this, at, args_str, error) my_local_virial_1 = at_local_virial_ptr if (store_contributions) call add_property(at, trim(calc_local_virial)//"_pot1", at_local_virial_ptr, overwrite=.true.) endif + if (len_trim(calc_local_dipole) > 0) then + call assign_property_pointer(at, trim(calc_local_dipole), at_local_dipole_ptr, error=error) + PASS_ERROR(error) + allocate(my_local_dipole_1(3, at%N)) + my_local_dipole_1 = at_local_dipole_ptr + if (store_contributions) call add_property(at, trim(calc_local_dipole)//"_pot1", at_local_dipole_ptr, overwrite=.true.) + endif call calc(this%pot2, at, args_str=trim(my_args_str)//" "//calc_args_pot2, error=error) @@ -141,6 +155,12 @@ recursive subroutine Potential_Sum_Calc(this, at, args_str, error) virial = my_virial_1 + virial call set_param_value(at, trim(calc_virial), virial) endif + if (len_trim(calc_dipole) > 0) then + call get_param_value(at, trim(calc_dipole), dipole) + if (store_contributions) call set_param_value(at, trim(calc_dipole)//"_pot2", dipole) + dipole = my_dipole_1 + dipole + call set_param_value(at, trim(calc_dipole), dipole) + endif if (len_trim(calc_local_energy) > 0) then if (store_contributions) call add_property(at, trim(calc_local_energy)//"_pot2", at_local_energy_ptr, overwrite=.true.) at_local_energy_ptr = my_local_e_1 + at_local_energy_ptr @@ -153,10 +173,16 @@ recursive subroutine Potential_Sum_Calc(this, at, args_str, error) if (store_contributions) call add_property(at, trim(calc_local_virial)//"_pot2", at_local_virial_ptr, overwrite=.true.) at_local_virial_ptr = my_local_virial_1 + at_local_virial_ptr end if + if (len_trim(calc_local_dipole) > 0) then + if (store_contributions) call add_property(at, trim(calc_local_dipole)//"_pot2", at_local_dipole_ptr, overwrite=.true.) + at_local_dipole_ptr = my_local_dipole_1 + at_local_dipole_ptr + end if + if (allocated(my_local_e_1)) deallocate(my_local_e_1) if (allocated(my_f_1)) deallocate(my_f_1) if (allocated(my_local_virial_1)) deallocate(my_local_virial_1) + if (allocated(my_local_dipole_1)) deallocate(my_local_dipole_1) end subroutine Potential_Sum_Calc diff --git a/src/Potentials/Potential_simple.F90 b/src/Potentials/Potential_simple.F90 index fe0bdb2aac..b43ac591be 100644 --- a/src/Potentials/Potential_simple.F90 +++ b/src/Potentials/Potential_simple.F90 @@ -358,7 +358,7 @@ recursive subroutine Potential_Simple_Calc(this, at, args_str, error) real(dp) :: energy, virial(3,3), deform(3,3), lat_save(3,3) real(dp), allocatable :: allpos_save(:,:) - real(dp), pointer :: at_force_ptr(:,:), at_local_energy_ptr(:), at_local_virial_ptr(:,:) + real(dp), pointer :: at_force_ptr(:,:), at_local_energy_ptr(:), at_local_virial_ptr(:,:) integer:: i,j,k,n, zero_loc(1) real(dp):: e_plus, e_minus, pos_save, r_scale, E_scale, cluster_box_buffer @@ -379,8 +379,10 @@ recursive subroutine Potential_Simple_Calc(this, at, args_str, error) integer, pointer :: cut_bonds_p(:,:), old_cut_bonds_p(:,:) integer :: i_inner, i_outer, n_non_term type(Atoms) :: cluster - character(len=STRING_LENGTH), target :: calc_force, calc_energy, calc_local_energy, calc_virial, calc_local_virial - logical :: do_calc_force, do_calc_energy, do_calc_local_energy, do_calc_virial, do_calc_local_virial + character(len=STRING_LENGTH), target :: calc_force, calc_energy, calc_local_energy, calc_virial, calc_local_virial, & + calc_dipole, calc_local_dipole + logical :: do_calc_force, do_calc_energy, do_calc_local_energy, do_calc_virial, do_calc_local_virial, & + do_calc_dipole, do_calc_local_dipole integer, pointer :: cluster_mark_p(:), at_prop_ptr_i(:), cluster_prop_ptr_i(:) integer, pointer :: old_cluster_mark_p(:) @@ -458,6 +460,10 @@ recursive subroutine Potential_Simple_Calc(this, at, args_str, error) help_string="If present, calculate virial and put it in field with this string as name") call param_register(params, 'local_virial', '', calc_local_virial, & help_string="If present, calculate local_virial and put it in field with this string as name") + call param_register(params, 'dipole', '', calc_dipole, & + help_string="If present, calculate dipole and put it in field with this string as name") + call param_register(params, 'local_dipole', '', calc_local_dipole, & + help_string="If present, calculate local_dipole and put it in field with this string as name") call param_register(params, "read_extra_param_list", '', read_extra_param_list, & help_string="if single_cluster=T and carve_cluster=T, extra params to copy back from cluster") call param_register(params, "read_extra_property_list", '', read_extra_property_list, & @@ -853,6 +859,9 @@ recursive subroutine Potential_Simple_Calc(this, at, args_str, error) do_calc_local_energy = len_trim(calc_local_energy) > 0 do_calc_virial = (len_trim(calc_virial) > 0) .and. .not. virial_using_fd do_calc_local_virial = len_trim(calc_local_virial) > 0 + do_calc_dipole = len_trim(calc_dipole) > 0 + do_calc_local_dipole = len_trim(calc_local_dipole) > 0 + call print("do_calc_force= "//do_calc_force//" calc_force="//trim(calc_force), PRINT_VERBOSE) call print("force_using_fd= "//force_using_fd, PRINT_VERBOSE) call print("virial_using_fd= "//virial_using_fd, PRINT_VERBOSE) @@ -860,6 +869,8 @@ recursive subroutine Potential_Simple_Calc(this, at, args_str, error) call print("do_calc_local_energy= "//do_calc_local_energy//" calc_local_energy="//trim(calc_local_energy), PRINT_VERBOSE) call print("do_calc_virial= "//do_calc_virial//" calc_virial="//trim(calc_virial), PRINT_VERBOSE) call print("do_calc_local_virial= "//do_calc_local_virial//" calc_local_virial="//trim(calc_local_virial), PRINT_VERBOSE) + call print("do_calc_dipole= "//do_calc_dipole//" calc_dipole="//trim(calc_dipole), PRINT_VERBOSE) + call print("do_calc_local_dipole= "//do_calc_local_dipole//" calc_local_dipole="//trim(calc_local_dipole), PRINT_VERBOSE) if(do_calc_virial .or. do_calc_energy .or. do_calc_force .or. do_calc_local_energy .or. do_calc_local_virial) then if(associated(this%ip)) then @@ -1286,6 +1297,12 @@ recursive subroutine Potential_Simple_Calc(this, at, args_str, error) end if + if( do_calc_dipole .or. do_calc_local_dipole ) then + if( associated(this%ip) ) then + call Calc(this%ip, at, args_str=args_str, error=error) + endif + endif + if (force_using_fd .and. len_trim(calc_force) > 0) then ! do forces by finite difference call print("Calculating force by finite differences with displacement="//force_fd_delta, PRINT_VERBOSE) From 3b16e7b04f3067fd529b3a44191fc4aaa66f3cc9 Mon Sep 17 00:00:00 2001 From: Albert Bartok-Partay Date: Sun, 5 Jul 2026 15:16:17 +0100 Subject: [PATCH 2/2] Update GAP version: dipole fitting in gap_fit --- src/GAP | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/GAP b/src/GAP index 480fc2d1c8..3ca04d2dda 160000 --- a/src/GAP +++ b/src/GAP @@ -1 +1 @@ -Subproject commit 480fc2d1c851d69eaa7932832f0c48c306bb5b93 +Subproject commit 3ca04d2dda03c0396b306a5f2129b1d13049d527