Skip to content
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.

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

At gaussian wt post mpi #3

Open
wants to merge 5 commits into
base: main
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from all commits
Commits
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
37 changes: 37 additions & 0 deletions descriptors.f95
Original file line number Diff line number Diff line change
Expand Up @@ -7140,6 +7140,11 @@ subroutine soap_calc(this,at,descriptor_out,do_descriptor,do_grad_descriptor,arg
complex(dp), allocatable, save :: sphericalycartesian_all_t(:,:), gradsphericalycartesian_all_t(:,:,:)
complex(dp) :: c_tmp(3)
integer :: max_n_neigh

character(len=STRING_LENGTH) :: atom_gaussian_weight_name
real(dp), dimension(:), pointer :: atom_gaussian_weight_pointer
logical :: has_atom_gaussian_weight_name

!$omp threadprivate(radial_fun, radial_coefficient, grad_radial_fun, grad_radial_coefficient)
!$omp threadprivate(sphericalycartesian_all_t, gradsphericalycartesian_all_t)
!$omp threadprivate(fourier_so3_r, fourier_so3_i)
Expand Down Expand Up @@ -7172,6 +7177,10 @@ subroutine soap_calc(this,at,descriptor_out,do_descriptor,do_grad_descriptor,arg
has_atom_mask_name = .false. ! allow atom mask column in the atom table
atom_mask_pointer => null() ! allow atom mask column in the atom table
xml_version = 1423143769 ! This is the version number where the 2l+1 normalisation of soap vectors was introduced

has_atom_gaussian_weight_name = .false.
atom_gaussian_weight_pointer => null()

if(present(args_str)) then
call initialise(params)

Expand All @@ -7182,6 +7191,10 @@ subroutine soap_calc(this,at,descriptor_out,do_descriptor,do_grad_descriptor,arg
call param_register(params, 'xml_version', '1423143769', xml_version, &
help_string="Version of GAP the XML potential file was created")

call param_register(params, 'atom_gaussian_weight_name', 'NONE', atom_gaussian_weight_name, &
has_value_target=has_atom_gaussian_weight_name, &
help_string="Array name from which to read prefactor for atom gaussians in the atomic neighbourhood density")

if (.not. param_read_line(params,args_str,ignore_unknown=.true.,task='soap_calc args_str')) then
RAISE_ERROR("soap_calc failed to parse args_str='"//trim(args_str)//"'", error)
endif
Expand All @@ -7196,6 +7209,17 @@ subroutine soap_calc(this,at,descriptor_out,do_descriptor,do_grad_descriptor,arg
atom_mask_pointer => null()
endif

if( has_atom_gaussian_weight_name ) then
if (.not. assign_pointer(at, trim(atom_gaussian_weight_name), atom_gaussian_weight_pointer)) then
RAISE_ERROR("soap_calc did not find "//trim(atom_gaussian_weight_name)//" property in the atoms object.", error)
endif
call print("got atom gaussian weight")
else
atom_mask_pointer => null() ! this would have been null anyway, so why to reassign it to null?
call print("DID NOT get atom gaussian weight")

endif

endif

if( this%cutoff_dexp > 0 ) then
Expand Down Expand Up @@ -7396,8 +7420,10 @@ subroutine soap_calc(this,at,descriptor_out,do_descriptor,do_grad_descriptor,arg
global_fourier_so3_i_array = 0.0_dp
endif ! this%global

! not sure if "atom_gaussian_weight_pointer has to be shared or if that's what messes things up?
!$omp parallel do schedule(dynamic) default(none) shared(this, at, descriptor_out, my_do_descriptor, my_do_grad_descriptor, d, i_desc, species_map, rs_index, do_two_l_plus_one) &
!$omp shared(global_grad_fourier_so3_r_array, global_grad_fourier_so3_i_array, norm_radial_decay) &
!$omp shared(atom_gaussian_weight_pointer) &
!$omp private(i, j, i_species, j_species, a, b, l, m, n, n_i, r_ij, u_ij, d_ij, shift_ij, i_pow, i_coeff, ia, jb, alpha, i_desc_i) &
!$omp private(c_tmp) &
!$omp private(t_g_r, t_g_i, t_f_r, t_f_i, t_g_f_rr, t_g_f_ii) &
Expand Down Expand Up @@ -7438,6 +7464,12 @@ subroutine soap_calc(this,at,descriptor_out,do_descriptor,do_grad_descriptor,arg
if( this%central_reference_all_species .or. this%species_Z(i_species) == at%Z(i) .or. this%species_Z(i_species) == 0 ) then
fourier_so3_r(0,a,i_species)%m(0) = this%central_weight * real(radial_coefficient(0,a) * SphericalYCartesian(0,0,(/0.0_dp, 0.0_dp, 0.0_dp/)), dp)
fourier_so3_i(0,a,i_species)%m(0) = this%central_weight * aimag(radial_coefficient(0,a) * SphericalYCartesian(0,0,(/0.0_dp, 0.0_dp, 0.0_dp/)))

if(associated(atom_gaussian_weight_pointer)) then
fourier_so3_r(0,a,i_species)%m(0) = fourier_so3_r(0,a,i_species)%m(0) * atom_gaussian_weight_pointer(i)
fourier_so3_i(0,a,i_species)%m(0) = fourier_so3_i(0,a,i_species)%m(0) * atom_gaussian_weight_pointer(i)
endif

else
fourier_so3_i(0,a,i_species)%m(0) = 0.0_dp
fourier_so3_r(0,a,i_species)%m(0) = 0.0_dp
Expand Down Expand Up @@ -7482,6 +7514,11 @@ subroutine soap_calc(this,at,descriptor_out,do_descriptor,do_grad_descriptor,arg
endif
f_cut = f_cut * radial_decay

if(associated(atom_gaussian_weight_pointer)) then
f_cut = f_cut * atom_gaussian_weight_pointer(j)
df_cut = df_cut * atom_gaussian_weight_pointer(j)
endif

do a = 1, this%n_max
arg_bess = 2.0_dp * this%alpha * r_ij * this%r_basis(a)
exp_p = exp( -this%alpha*( r_ij + this%r_basis(a) )**2 )
Expand Down
19 changes: 17 additions & 2 deletions gap_fit_module.f95
Original file line number Diff line number Diff line change
Expand Up @@ -116,6 +116,8 @@ module gap_fit_module
logical :: sparse_use_actual_gpcov
logical :: has_template_file, has_e0, has_local_property0, has_e0_offset, has_linear_system_dump_file

character(len=STRING_LENGTH) :: descriptor_args_str

endtype gap_fit

private
Expand Down Expand Up @@ -180,6 +182,12 @@ subroutine gap_fit_parse_command_line(this)
integer :: rnd_seed
integer, pointer :: mpi_blocksize

! this was the bug, I think
! character(len=STRING_LENGTH) :: descriptor_args_str
character(len=STRING_LENGTH), pointer :: descriptor_args_str
logical :: has_descriptor_args_str


config_file => this%config_file
at_file => this%at_file
e0_str => this%e0_str
Expand Down Expand Up @@ -208,6 +216,8 @@ subroutine gap_fit_parse_command_line(this)
gp_file => this%gp_file
template_file => this%template_file
sparsify_only_no_fit => this%sparsify_only_no_fit
descriptor_args_str => this%descriptor_args_str

condition_number_norm => this%condition_number_norm
linear_system_dump_file => this%linear_system_dump_file
mpi_blocksize => this%mpi_blocksize
Expand Down Expand Up @@ -331,9 +341,12 @@ subroutine gap_fit_parse_command_line(this)
call param_register(params, 'template_file', 'template.xyz', template_file, has_value_target=this%has_template_file, &
help_string="Template XYZ file for initialising object")

call param_register(params, 'descriptor_args_str', '', descriptor_args_str, &
help_string="Arguments string for descriptor")

call param_register(params, 'sparsify_only_no_fit', 'F', sparsify_only_no_fit, &
help_string="If true, sparsification is done, but no fitting. print the sparse index by adding print_sparse_index=file.dat to the descriptor string.")

call param_register(params, 'condition_number_norm', ' ', condition_number_norm, &
help_string="Norm for condition number of matrix A; O: 1-norm, I: inf-norm, <space>: skip calculation (default)")

Expand Down Expand Up @@ -1217,7 +1230,9 @@ subroutine fit_data_from_xyz(this,error)
do i_coordinate = 1, this%n_coordinate

call calc(this%my_descriptor(i_coordinate),this%at(n_con),my_descriptor_data, &
do_descriptor=.true.,do_grad_descriptor=has_force .or. has_virial)
do_descriptor=.true.,do_grad_descriptor=has_force .or. has_virial, args_str=this%descriptor_args_str)
! old: check this if doesn't work
! do_descriptor=.true.,do_grad_descriptor=has_force .or. has_virial)

allocate(xloc(size(my_descriptor_data%x)))
n_descriptors = n_descriptors + size(my_descriptor_data%x)
Expand Down