Visualisation

.gitignore

@@ -28,3 +28,5 @@ fort.*
 *.e
 *.log
 *.agr
+*.pdf
+*.eps

example/wrapper/run_BaTiO3.py

@@ -27,7 +27,7 @@
 
 # read the host structure from a POSCAR file
 print("Reading host")
-host = read("../example_files/POSCAR_host_diamond")
+host = read("../example_files/POSCAR_host_BaTiO3")
 host.calc = calculator
 print("host energy: ", host.get_potential_energy())
 
@@ -137,11 +137,11 @@
 # generate structures
 num_structures_old = 0
 optimise_structure = False
-for iter in range(3):
+for iter in range(1):
     print(f"Iteration {iter}")
     print("Generating...")
     # this is the main function to generate structures
-    generator.generate(num_structures=1, stoichiometry=stoich_dict, seed=4+iter, verbose=1, method_probab={"void":0.001, "walk":0.0, "min":1.0})
+    generator.generate(num_structures=1, stoichiometry=stoich_dict, seed=0+iter, verbose=1, method_probab={"void":0.001, "walk":0.0, "min":1.0})
     print("Generated")
 
     print("Getting structures")

src/fortran/lib/mod_evaluator.f90

@@ -221,7 +221,6 @@ function evaluate_point( gvector_container, &
        viability_2body = 0.5_real12
     else
        viability_2body = viability_2body / real( num_2body, real12 )
-      !  viability_2body = viability_2body ** (1._real12 / num_2body)
     end if
 
 
@@ -275,15 +274,17 @@ function evaluate_point( gvector_container, &
     ! Normalise the viability map
     if(num_3body.eq.0)then
        viability_3body = gvector_container%viability_3body_default
-      !  viability_3body = 1.5_real12
     else
-       viability_3body = viability_3body ** (1._real12 / real(num_3body,real12))
+       viability_3body = viability_3body ** ( &
+            1._real12 / real(num_3body,real12) &
+       )
     end if
     if(num_4body.eq.0)then
        viability_4body = gvector_container%viability_4body_default
-      !  viability_4body = 1.5_real12
     else
-       viability_4body = viability_4body ** (1._real12 / real(num_4body,real12))
+       viability_4body = viability_4body ** ( &
+            1._real12 / real(num_4body,real12) &
+       )
     end if
 
     ! Combine the 2-, 3- and 4-body maps
@@ -326,11 +327,10 @@ function evaluate_3body_contributions( gvector_container, &
 
 
     output = 1._real12
-    num_3body_local = 0
+    num_3body_local = sum(basis%spec(current_idx(1):)%num) - current_idx(2)
     species_loop: do js = current_idx(1), basis%nspec, 1
        atom_loop: do ja = 1, basis%spec(js)%num
-          if(all([js,ja].eq.current_idx))cycle
-          num_3body_local = num_3body_local + 1
+          if(js.eq.current_idx(1) .and. ja.le.current_idx(2))cycle
           associate( position_store => [ basis%spec(js)%atom(ja,1:3) ] )
              bin = gvector_container%get_bin( &
                   get_angle( position_2, &
@@ -345,15 +345,17 @@ function evaluate_3body_contributions( gvector_container, &
                 write(0,*) "Error: bin = 0, IF NOT TRIGGERED, WE CAN REMOVE THIS IF"
                 stop 1
              end if
-             output = output * gvector_container%total%df_3body(bin,species)
+             output = output * &
+                  gvector_container%total%df_3body( &
+                       bin, &
+                       gvector_container%host_system%element_map(species) &
+                  ) ** ( 1._real12 / real( num_3body_local, real12 ) )
           end associate
        end do atom_loop
     end do species_loop
     if(num_3body_local.eq.0)then
        output = 1._real12
        num_3body = num_3body - 1
-    else
-       output = output ** (1._real12 / real( num_3body_local, real12 ) )
     end if
 
   end function evaluate_3body_contributions
@@ -410,9 +412,10 @@ function evaluate_4body_contributions( gvector_container, &
                 stop 1
              end if
              output = output * &
-                  gvector_container%total%df_4body(bin,species) ** ( &
-                       1._real12 / real( num_4body_local, real12 ) &
-                  )
+                  gvector_container%total%df_4body( &
+                       bin, &
+                       gvector_container%host_system%element_map(species) &
+                  ) ** ( 1._real12 / real( num_4body_local, real12 ) )
           end associate
        end do atom_loop
     end do species_loop

src/fortran/lib/mod_evolver.f90

@@ -75,11 +75,16 @@ module evolver
      type(basis_type) :: basis
      !! Host structure.
      integer, dimension(:,:), allocatable :: pair_index
+     !! Index for the 2-body distribution function.
+     integer, dimension(:), allocatable :: element_map
+     !! Mapping of host elements to distribution function elements.
    contains
      procedure, pass(this) :: calculate_interface_energy
      !! Calculate the interface formation energy of the host.
      procedure, pass(this) :: set => set_host
      !! Set the host structure for the distribution functions.
+     procedure, pass(this) :: set_element_map => set_host_element_map
+     !! Set the mapping of host elements to distribution function elements.
   end type gvector_host_type
 
   type :: gvector_container_type
@@ -97,6 +102,8 @@ module evolver
      !! Boolean whether to weight the distribution functions by the energy
      !! above the hull. If false, the formation energy from the element
      !! reference energies is used.
+     integer, dimension(:), allocatable :: host_to_df_species_map
+     !! Mapping of host species to distribution function species.
      real(real12) :: &
           viability_3body_default = 0.1_real12, &
           viability_4body_default = 0.1_real12
@@ -222,6 +229,8 @@ module evolver
      !! Write the learned 4-body distribution function to a file.
      procedure, pass(this) :: get_pair_index
      !! Return the index for bond_info given two elements.
+     procedure, pass(this) :: get_element_index
+     !! Return the index for element_info given one element.
      procedure, pass(this) :: get_bin
      !! Return the bin index for a given distance.
   end type gvector_container_type
@@ -431,6 +440,12 @@ subroutine set_host(this, host)
    if(allocated(this%df_3body)) deallocate(this%df_3body)
    if(allocated(this%df_4body)) deallocate(this%df_4body)
 
+   !! Run set_element_map if total dfs have already been calculated
+   !! else, this will be run in the create() procedure
+   if(allocated(this%total%df_2body))then
+      call this%set_element_map(this%element_info)
+   end if
+
   end subroutine set_host
 !###############################################################################
 
@@ -551,6 +566,8 @@ subroutine create( &
     call this%set_bond_info()
     call this%evolve()
     if(deallocate_systems_) call this%deallocate_systems()
+    if(this%host_system%defined) &
+         call this%host_system%set_element_map(this%element_info)
 
   end subroutine create
 !###############################################################################
@@ -673,6 +690,8 @@ subroutine update( &
 
     call this%evolve()
     if(deallocate_systems_) call this%deallocate_systems()
+    if(this%host_system%defined) &
+         call this%host_system%set_element_map(this%element_info)
 
   end subroutine update
 !###############################################################################
@@ -1798,6 +1817,61 @@ end function get_pair_index
 !###############################################################################
 
 
+!###############################################################################
+  pure function get_element_index(this, species) result(idx)
+    !! Get the index of an element in the container.
+    implicit none
+
+    ! Arguments
+    class(gvector_container_type), intent(in) :: this
+    !! Parent of the procedure. Instance of distribution functions container.
+    character(len=3), intent(in) :: species
+    !! Element name.
+    integer :: idx
+    !! Index of the element in the element_info array.
+
+    ! Local variables
+    integer :: is, js
+    !! Index of the elements in the element_info array.
+
+    idx = findloc([ this%element_info(:)%name ], species, dim=1)
+
+  end function get_element_index
+!###############################################################################
+
+
+!###############################################################################
+  subroutine set_host_element_map(this, element_info)
+    !! Set the host element map for the container.
+    implicit none
+
+    ! Arguments
+    class(gvector_host_type), intent(inout) :: this
+    !! Parent of the procedure. Instance of distribution functions container.
+    type(element_type), dimension(:), intent(in) :: element_info
+    !! Element information.
+
+    ! Local variables
+    integer :: is, js
+    !! Index of the elements in the element_info array.
+
+    if(.not.this%defined)then
+       call stop_program( "Host not defined" )
+       return
+    end if
+    if(allocated(this%element_map)) deallocate(this%element_map)
+    allocate(this%element_map(this%basis%nspec))
+    do is = 1, this%basis%nspec
+       this%element_map(is) = findloc(&
+            [ element_info(:)%name ], &
+            this%basis%spec(is)%name, dim=1 &
+       )
+    end do
+
+  end subroutine set_host_element_map
+!###############################################################################
+
+
 !###############################################################################
   pure function get_bin(this, value, dim) result(bin)
     !! Get the bin index for a value in a dimension.

src/fortran/lib/mod_generator.f90

@@ -154,7 +154,7 @@ subroutine set_host(this, host)
     ! Arguments
     class(raffle_generator_type), intent(inout) :: this
     !! Instance of the raffle generator.
-    type(basis_type), intent(in) :: host
+    class(basis_type), intent(in) :: host
     !! Basis of the host structure.
 
     ! Local variables

src/fortran/lib/mod_misc.f90

@@ -1069,6 +1069,11 @@ end function to_lower
 !###############################################################################
   function strip_null(buffer) result(stripped)
     !! Strip null characters from a string.
+    !!
+    !! This is meant for handling strings passed from Python, which gain
+    !! null characters at the end. The procedure finds the first null
+    !! character and truncates the string at that point.
+    !! Null characters are represented by ASCII code 0.
     implicit none
 
     ! Arguments

src/fortran/lib/mod_misc_linalg.f90

@@ -107,10 +107,15 @@ pure function get_angle_from_vectors(vector1, vector2) result(angle)
     real(real12) :: angle
     !! Output angle.
 
-    angle = acos( dot_product(vector1,vector2)/&
-         ( modu(vector1) * modu(vector2) ))
-    if (isnan(angle)) angle = 0._real12
-
+    angle =  dot_product(vector1,vector2) / &
+         ( modu(vector1) * modu(vector2) )
+    if(angle .ge. 1._real12)then
+       angle = 0._real12
+    elseif(angle .le. -1._real12)then
+       angle = pi
+    else
+       angle = acos(angle)
+    end if
   end function get_angle_from_vectors
 !###############################################################################
 
@@ -128,9 +133,15 @@ pure function get_angle_from_points(point1, point2, point3) result(angle)
     real(real12) :: angle
     !! Output angle.
 
-    angle = acos( ( dot_product( point1 - point2, point3 - point2 ) ) / &
-         ( modu( point1 - point2 ) * modu( point3 - point2 ) ) )
-    if(isnan(angle)) angle = 0._real12
+    angle = dot_product( point1 - point2, point3 - point2 ) / &
+         ( modu( point1 - point2 ) * modu( point3 - point2 ) )
+    if(angle .ge. 1._real12)then
+       angle = 0._real12
+    elseif(angle .le. -1._real12)then
+       angle = pi
+    else
+       angle = acos(angle)
+    end if
   end function get_angle_from_points
 !###############################################################################
 

test/CMakeLists.txt

@@ -9,7 +9,8 @@ foreach(execid
     extended_geom
     atom_adder
     evolver
-    evaluator
+    evaluator_C
+    evaluator_BTO
     generator
   )
   add_executable(test_${execid} test_${execid}.f90)