Add capability for electrostatics from ML charges

.gitignore

@@ -2,3 +2,5 @@ bin/
 build/
 include/
 lib/
+trajectory_out.xyz
+thermo.log

Makefile

@@ -23,7 +23,7 @@ PROGRAMS := turbogap
 
 SRC := splines.f90 types.f90 neighbors.f90 gap.f90 vdw.f90		\
 	local_properties.f90 exp_utils.f90  xyz.f90 md.f90 mc.f90 read_files.f90	\
-	gap_interface.f90 mpi.f90 exp_interface.f90
+	gap_interface.f90 mpi.f90 exp_interface.f90 electrostatics.f90
 
 SRC_TP_BT := resamplekin.f90
 SRC_ST := soap_turbo_functions.f90 soap_turbo_radial.f90 soap_turbo_angular.f90 \

src/gap_interface.f90

@@ -277,7 +277,8 @@ subroutine get_gap_soap(n_total_sites, n_sites0, n_neigh0, neighbors_list0, n_sp
                & local_property_models(i4)%delta,&
                & local_property_models(i4)%zeta, local_properties,&
                & do_derivatives, soap_cart_der, n_neigh,&
-               & local_properties_cart_der )
+               & local_properties_cart_der, &
+               & local_property_models(i4)%zero_trunc)
 
           ! call local_property_predict( soap, Qs, alphas, V0, delta, zeta, &
           !      local_property, do_derivatives, soap_cart_der, n_neigh_out, &
@@ -766,15 +767,15 @@ subroutine get_local_properties( soap, Qs, alphas, V0, delta, zeta, &
                               local_property0, do_derivatives, soap_cart_der, n_neigh_out, &
                               local_property_cart_der0, n_atom_pairs,&
                               & in_to_out_pairs, n_all_sites,&
-                              & in_to_out_site, n_sites_out )
+                              & in_to_out_site, n_sites_out, do_zero_floor )
 
   real*8, allocatable, intent(in) :: soap(:,:), soap_cart_der(:,:,:)
   real*8, intent(in) ::  Qs(:,:), alphas(:), zeta, delta, V0
   real*8, intent(inout) :: local_property0(:), local_property_cart_der0(:,:)
   real*8, allocatable :: local_property(:), local_property_cart_der(:,:)
   integer, intent(in) :: n_atom_pairs, in_to_out_pairs(:), n_all_sites,&
        & in_to_out_site(:), n_neigh_out(:), n_sites_out
-  logical, intent(in) :: do_derivatives
+  logical, intent(in) :: do_derivatives, do_zero_floor
 
   allocate( local_property( 1:n_sites_out ) )
   local_property = 0.d0
@@ -785,7 +786,7 @@ subroutine get_local_properties( soap, Qs, alphas, V0, delta, zeta, &
 
   call local_property_predict( soap, Qs, alphas, V0, delta, zeta, &
        local_property, do_derivatives, soap_cart_der, n_neigh_out, &
-       local_property_cart_der )
+       local_property_cart_der, do_zero_floor)
   do i = 1, n_sites_out
      i2 = in_to_out_site(i)
      local_property0(i2) = local_property(i)

src/local_properties.f90

@@ -176,15 +176,16 @@ module local_prop
     ! end subroutine get_local_property_details
 
 
-  subroutine local_property_predict( soap, Qs, alphas, V0, delta, zeta, V, &
-                                do_derivatives, soap_cart_der, n_neigh, V_der )
+  subroutine local_property_predict(soap, Qs, alphas, V0, delta, zeta, V, &
+                                    do_derivatives, soap_cart_der, n_neigh, V_der, &
+                                    do_zero_floor)
 
     implicit none
 
 !   Input variables
     real*8, intent(in) :: soap(:,:), Qs(:,:), alphas(:), V0, delta, zeta, soap_cart_der(:,:,:)
     integer, intent(in) :: n_neigh(:)
-    logical, intent(in) :: do_derivatives
+    logical, intent(in) :: do_derivatives, do_zero_floor
 !   Output variables
     real*8, intent(out) :: V(:), V_der(:,:)
 !   Internal variables
@@ -249,12 +250,14 @@ subroutine local_property_predict( soap, Qs, alphas, V0, delta, zeta, V, &
     end if
     V = V + V0
 
-!   Make sure all V are >= 0
-    do i = 1, size(V)
-      if( V(i) < 0.d0 )then
-        V(i) = 0.d0
-      end if
-    end do
+    !   Make sure all V are >= 0, if requested
+    if (do_zero_floor) then
+      do i = 1, size(V)
+        if( V(i) < 0.d0 )then
+          V(i) = 0.d0
+        end if
+      end do
+    end if
 
     if( do_derivatives)then
       if( n_sites > 0 )then

src/read_files.f90

@@ -975,6 +975,9 @@ subroutine read_input_file(n_species, mode, params, rank)
       else if(keyword=='print_vdw_forces')then
         backspace(10)
         read(10, *, iostat=iostatus) cjunk, cjunk, params%print_vdw_forces
+      else if(keyword=='print_estat_forces')then
+        backspace(10)
+        read(10, *, iostat=iostatus) cjunk, cjunk, params%print_estat_forces
       else if(keyword=='exp_similarity_type')then
         backspace(10)
         read(10, *, iostat=iostatus) cjunk, cjunk, params%exp_similarity_type
@@ -1271,9 +1274,11 @@ subroutine read_input_file(n_species, mode, params, rank)
          backspace(10)
          read(10, *, iostat=iostatus) cjunk, cjunk, &
               (params%exp_data(nw)%n_samples, nw=1, params%n_exp)
-
-      else if (keyword(i2-4:i2) == "range" .or.  keyword(i2-8:i2) ==&
-           & "file_data" .or. keyword(i2-8:i2) == "n_samples"  )then
+      !TODO WARNING this sort of thing is fragile and is causing unpredictable failures!
+      ! fixed it for now, but... there has to be a better way.
+      ! Update: the index() function would be cleaner, if a tad less efficient
+      else if (keyword(max(1,i2-4):i2) == "range" .or.  keyword(max(1, i2-8):i2) ==&
+           & "file_data" .or. keyword(max(1,i2-8):i2) == "n_samples"  )then
          backspace(10)
          ! Check if experimental range or data files are specified
          do nw = 1, params%n_exp
@@ -1654,6 +1659,29 @@ subroutine read_input_file(n_species, mode, params, rank)
       else if( keyword == "print_vdw_forces" )then
         backspace(10)
         read(10, *, iostat=iostatus) cjunk, cjunk, params%print_vdw_forces
+      else if( keyword == "estat_method" )then
+        backspace(10)
+        read(10, *, iostat=iostatus) cjunk, cjunk, params%estat_method
+        params%estat_method = trim(params%estat_method)
+        call upper_to_lower_case(params%estat_method)
+        if (params%estat_method /= "direct" .and. params%estat_method /= "dsf" &
+            .and. params%estat_method /= "none") then
+          write(*,*) "ERROR: electrostatic method not implemented: ", params%estat_method
+          write(*,*) "       Currently implemented methods are: direct"
+          stop
+        end if
+      else if( keyword == "estat_dsf_alpha" )then
+        backspace(10)
+        read(10, *, iostat=iostatus) cjunk, cjunk, params%estat_dsf_alpha
+      else if( keyword == "estat_rcut" )then
+        backspace(10)
+        read(10, *, iostat=iostatus) cjunk, cjunk, params%estat_rcut
+      else if( keyword == "estat_rcut_inner" )then
+        backspace(10)
+        read(10, *, iostat=iostatus) cjunk, cjunk, params%estat_rcut_inner
+      else if( keyword == "estat_inner_width" )then
+        backspace(10)
+        read(10, *, iostat=iostatus) cjunk, cjunk, params%estat_inner_width
       else if( keyword == "core_pot_cutoff" )then
         backspace(10)
         read(10, *, iostat=iostatus) cjunk, cjunk, params%core_pot_cutoff
@@ -1739,6 +1767,14 @@ subroutine read_input_file(n_species, mode, params, rank)
       end if
     end if
 
+    ! Do the same for electrostatics
+    if (params%estat_method == "none") then
+      params%estat_rcut = 0.0_dp
+    else if (params%estat_method == "dsf" .and. params%estat_dsf_alpha <= 0) then
+      write(*,*) "ERROR: Must specify an alpha (>0) for DSF electrostatics."
+      stop
+    end if
+
 !   Nested sampling checks
     if( params%do_nested_sampling )then
       if( params%thermostat /= "none" )then
@@ -2340,6 +2376,10 @@ subroutine read_gap_hypers(file_gap, &
                     soap_turbo_hypers(n_soap_turbo)%core_electron_be_index=nw
                  end if
 
+                 if (trim(soap_turbo_hypers(n_soap_turbo)%local_property_models(nw)%label) &
+                        == "atomic_charge") then
+                    soap_turbo_hypers(n_soap_turbo)%has_charges = .true.
+                 end if
               end do
 
             else if( keyword == "local_property_qs" )then
@@ -2799,6 +2839,189 @@ subroutine upper_to_lower_case(string)
 !**************************************************************************
 
 
+  subroutine get_irreducible_local_properties(params, n_local_properties_tot, n_soap_turbo, soap_turbo_hypers, &
+       local_property_labels, local_property_labels_temp, local_property_labels_temp2, local_property_indexes, &
+       valid_vdw, vdw_lp_index, valid_estat_charges, charge_lp_index, core_be_lp_index, valid_xps, xps_idx )
+    implicit none
+    type( input_parameters ), intent(inout) :: params
+    integer, intent(in) :: n_soap_turbo
+    integer, intent( inout ) :: n_local_properties_tot
+    type( soap_turbo ), allocatable, intent(inout) :: soap_turbo_hypers(:)
+    character*1024, allocatable, intent(inout) ::  local_property_labels(:), local_property_labels_temp(:), &
+         local_property_labels_temp2(:)
+    integer, allocatable, intent(inout) :: local_property_indexes(:)
+    integer, intent(inout) :: vdw_lp_index, core_be_lp_index, charge_lp_index, xps_idx
+    logical, intent(inout) :: valid_vdw, valid_xps, valid_estat_charges
+    logical :: label_in_list = .false.
+    integer :: i, j, i2, j2, k, k2, nprop
+
+    n_local_properties_tot = 0
+    i2 = 1 ! using this as a counter for the labels
+    do j = 1, n_soap_turbo
+       if( soap_turbo_hypers(j)%has_local_properties )then
+          ! This property has the labels of the quantities to
+          ! compute. We must specify the number of local properties, for the sake of coding simplicity
+
+          n_local_properties_tot = n_local_properties_tot + soap_turbo_hypers(j)%n_local_properties
+
+          if(.not. allocated(local_property_labels))then
+             allocate(local_property_labels(1:n_local_properties_tot))
+             do i = 1, n_local_properties_tot
+                local_property_labels(i) = soap_turbo_hypers(j)%local_property_models(i)%label
+                write(*,*)' Local property found                  |'
+                write(*,'(A,1X,I8,1X,A,1X,A)')' Descriptor ', j,&
+                     & trim(soap_turbo_hypers(j)&
+                     &%local_property_models(i)%label),  ' |'
+             end do
+          else
+             ! Allocate temporary array which is of the size before 
+             allocate( local_property_labels_temp( 1:n_local_properties_tot - soap_turbo_hypers(j)%n_local_properties ))
+             local_property_labels_temp = local_property_labels
+             deallocate(local_property_labels)
+             allocate(local_property_labels(1:n_local_properties_tot))
+
+             nprop = soap_turbo_hypers(j)%n_local_properties
+             do i = 1, n_local_properties_tot - nprop
+                local_property_labels(i) = local_property_labels_temp(i)
+             end do
+
+             deallocate(local_property_labels_temp)
+
+             do i = 1, nprop
+                local_property_labels(i + n_local_properties_tot -&
+                     & nprop) = soap_turbo_hypers(j)&
+                     &%local_property_models(i)%label
+                write(*,*)' Local property found                  |'
+                write(*,'(A,1X,I8,1X,A,1X,A)')' Descriptor ', j,&
+                     & trim(soap_turbo_hypers(j)&
+                     &%local_property_models(i)%label),  ' |'
+
+             end do
+          end if
+       end if
+    end do
+
+    ! by this point, local_property_labels( 1:n_local_properties_tot ) has labels of all local properties 
+
+
+    ! Now we create an irreducible list of the labels
+    i2=0
+    if (n_local_properties_tot > 0)then
+       allocate( local_property_labels_temp( 1:1 ))
+       local_property_labels_temp(1) = local_property_labels(1)
+       i2 = 1
+       if (n_local_properties_tot > 1)then
+          do i = 2, n_local_properties_tot
+             label_in_list = .false.
+             ! Iterate through irreducible list to see if there is a mismatch
+             do j = 1, size( local_property_labels_temp, 1 )
+                if (trim( local_property_labels_temp(j) ) == trim( local_property_labels(i) )) label_in_list = .true.
+             end do
+             if (.not. label_in_list) then
+                i2 = i2 + 1
+                allocate(local_property_labels_temp2(1:i2))
+                local_property_labels_temp2(1:i2-1) = local_property_labels_temp(1:i2-1)
+                local_property_labels_temp2(i2)     = local_property_labels(i)
+                deallocate(local_property_labels_temp)
+                allocate(local_property_labels_temp(1:i2))
+                local_property_labels_temp(1:i2) = local_property_labels_temp2(1:i2)
+                deallocate(local_property_labels_temp2)
+             end if
+          end do
+       end if
+
+       params%n_local_properties = i2
+
+       ! by this point, local_property_labels( 1:n_local_properties_tot ) has labels of all local properties
+       !                local_property_labels_temp( 1:params%n_local_properties ) has irreducible labels of local properties 
+
+       ! Now we can have an array which has a soap turbo index as an input and it can give us the corresponding label
+       allocate(local_property_indexes(1:n_local_properties_tot))
+       i2 = 1
+       do i = 1, params%n_local_properties
+          do j = 1, n_local_properties_tot
+             if ( trim(local_property_labels(j)) == trim( local_property_labels_temp(i) ) )then
+
+                local_property_indexes(j) = i
+
+                ! WARNING -- this seems to be duplicating functionality found above
+                ! see the line:
+                !   else if( keyword == "local_property_labels" )then
+                ! also perhaps in turbogap.f90, just below the get_gap_soap() call
+                if ( trim(local_property_labels(j)) == "hirshfeld_v" )then
+                   vdw_lp_index = i
+                   valid_vdw = .true.
+                   do k2 = 1, n_soap_turbo
+                      do k = 1, soap_turbo_hypers(k2)%n_local_properties
+                         if (trim(soap_turbo_hypers(k2)%local_property_models(k)%label) == "hirshfeld_v")then
+                            if( params%do_derivatives .or. params%do_forces)then
+                               soap_turbo_hypers(k2)%local_property_models(k)%do_derivatives = .true.
+                            else
+                               soap_turbo_hypers(k2)%local_property_models(k)%do_derivatives = .false.
+                            end if
+                            ! The previous default behaviour is to truncate negative
+                            ! volumes to zero, so we keep it that way here for
+                            ! backwards compatibility
+                            soap_turbo_hypers(k2)%local_property_models(k)%zero_trunc = .true.
+                         end if
+                      end do
+                   end do
+                end if
+
+                if (trim(local_property_labels(j)) == "atomic_charge") then
+                   charge_lp_index = i
+                   valid_estat_charges = .true.
+                   do k2 = 1, n_soap_turbo
+                      do k = 1, soap_turbo_hypers(k2)%n_local_properties
+                         if (trim(soap_turbo_hypers(k2)%local_property_models(k)%label) == "atomic_charge")then
+                            soap_turbo_hypers(k2)%has_charges = .true.
+                            if( params%do_derivatives .or. params%do_forces)then
+                               soap_turbo_hypers(k2)%local_property_models(k)%do_derivatives = .true.
+                            else
+                               soap_turbo_hypers(k2)%local_property_models(k)%do_derivatives = .false.
+                            end if
+                            ! This is important -- no truncating the charges; it doesn't make sense here!
+                            soap_turbo_hypers(k2)%local_property_models(k)%zero_trunc = .false.
+                         end if
+                      end do
+                   end do
+                end if
+
+                if ( trim(local_property_labels(j)) == "core_electron_be" )then
+                   core_be_lp_index = i
+
+                   ! Check if there is experimental data for one to do xps fitting
+                   do i2 = 1, params%n_exp
+                      if(( trim(params%exp_data(i2)%label) == "xps" .and.  &
+                           .not. ( trim(params%exp_data(i2)%file_data) == "none" )))then
+                         valid_xps = .true.
+                         xps_idx = i2
+                         do k2 = 1, n_soap_turbo
+                            do k = 1, soap_turbo_hypers(k2)%n_local_properties
+                               if (trim(soap_turbo_hypers(k2)%local_property_models(k)%label) == "core_electron_be")then
+                                  soap_turbo_hypers(k2)%local_property_models(k)%do_derivatives = .false.
+                                  if( params%exp_forces .and. params%do_derivatives)then
+                                     soap_turbo_hypers(k2)%local_property_models(k)%do_derivatives = .true.
+                                  end if
+                               end if
+                            end do
+                         end do
+                      end if
+                   end do
+                end if
+
+             end if
+          end do
+       end do
+
+       deallocate(local_property_labels)
+       allocate(local_property_labels(1:size(local_property_labels_temp,1)))
+       local_property_labels = local_property_labels_temp
+       deallocate(local_property_labels_temp)
+    end if
+
+  end subroutine get_irreducible_local_properties
+
 
 
 end module