Skip to content

Dipole fitting and prediction #727

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
albapa wants to merge 2 commits into
base: public
Choose a base branch
from
Open

Dipole fitting and prediction #727

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
e22cb24
Dipole and local dipole prediction added.
albapa Jul 5, 2026
3b16e7b
Update GAP version: dipole fitting in gap_fit
albapa Jul 5, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion src/GAP
Open in desktop
Submodule GAP updated 2 files
+19 −22 descriptors.F90
+143 −42 gap_fit_module.F90
65 changes: 56 additions & 9 deletions src/Potentials/IPModel_GAP.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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, &
Expand All @@ -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)//"'")
Expand Down Expand Up @@ -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')
Expand All @@ -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)
Expand All @@ -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
Expand All @@ -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)

Expand Down Expand Up @@ -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
Expand Down
36 changes: 31 additions & 5 deletions src/Potentials/Potential_Sum_routines.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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")
Expand All @@ -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)
Expand All @@ -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)
Expand All @@ -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
Expand All @@ -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

Expand Down
23 changes: 20 additions & 3 deletions src/Potentials/Potential_simple.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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(:)
Expand Down Expand Up @@ -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, &
Expand Down Expand Up @@ -853,13 +859,18 @@ 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)
call print("do_calc_energy= "//do_calc_energy//" calc_energy="//trim(calc_energy), PRINT_VERBOSE)
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
Expand Down Expand Up @@ -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)
Expand Down
Loading