From 66645c8dfc79108565e5b43931888d8655a4d456 Mon Sep 17 00:00:00 2001
From: Alex Selimov <alex@alexselimov.xyz>
Date: Fri, 18 Jun 2021 10:04:54 -0400
Subject: [PATCH] Added opt_bubble

---
 src/Makefile.dep   |   5 ++
 src/caller.f90     |   3 ++
 src/opt_bubble.f90 | 116 +++++++++++++++++++++++++++++++++++++++++++++
 src/opt_delete.f90 |   2 +-
 src/parameters.f90 |   1 +
 5 files changed, 126 insertions(+), 1 deletion(-)
 create mode 100644 src/opt_bubble.f90

diff --git a/src/Makefile.dep b/src/Makefile.dep
index d4f2675..af115e0 100644
--- a/src/Makefile.dep
+++ b/src/Makefile.dep
@@ -122,6 +122,11 @@ obj/neighbors.o :  \
 	obj/parameters.o \
 	obj/subroutines.o
 
+obj/opt_bubble.o :  \
+	obj/box.o \
+	obj/elements.o \
+	obj/parameters.o
+
 obj/opt_deform.o :  \
 	obj/box.o \
 	obj/elements.o \
diff --git a/src/caller.f90 b/src/caller.f90
index 18e91cf..dc99fbd 100644
--- a/src/caller.f90
+++ b/src/caller.f90
@@ -15,6 +15,7 @@ module caller
     use opt_delete
     use opt_redef_box
     use opt_slip_plane
+    use opt_bubble
     use box
 
 
@@ -87,6 +88,8 @@ module caller
             call redef_box(arg_pos)
         case('-slip_plane')
             call run_slip_plane(arg_pos)
+        case('-bubble')
+            call bubble(arg_pos)
         case default
             print *, 'Option ', trim(adjustl(option)), ' is not currently accepted.'
             stop 3 
diff --git a/src/opt_bubble.f90 b/src/opt_bubble.f90
new file mode 100644
index 0000000..5eaf2f5
--- /dev/null
+++ b/src/opt_bubble.f90
@@ -0,0 +1,116 @@
+module opt_bubble
+    !This module contains the bubble option which can be used to create voids with specific pressures of certain atoms
+
+    use parameters
+    use elements 
+    use box
+    use opt_group
+
+    implicit none
+
+    real(kind=dp), private :: br, bt, bp, c(3)
+    character(len=2), private :: species
+
+    public 
+    contains
+    subroutine bubble(arg_pos)
+        integer, intent(inout) :: arg_pos
+
+        integer :: new_type, n, j, i
+        real(kind = dp)  :: p(3), rand, factor, per, vol
+
+        print *, '------------------------------------------------------------'
+        print *, 'Option Bubble'
+        print *, '------------------------------------------------------------'
+        !First we parse the bubble command
+        call parse_bubble(arg_pos)
+
+        !Now we use the existing group code to delete a sphere which represents the bubble
+        centroid=c
+        radius = br
+        type='both'
+        gshape='sphere'
+        group_nodes = .true.
+
+        call get_group
+        call refine_group
+        call get_group
+        call delete_group
+
+        !Now we create a new atom type with the desired species
+        call add_atom_type(species, new_type)
+
+        !Now we calculate the number of atoms we need to add for the desired pressure
+        print *, "Creating a bubble with pressure", bp, " at temperature ", bt, " with radius ", br
+        
+        factor = 1.0e24/6.02214e23
+        if (bp <= 10.0) then
+            per=factor*(3.29674113+4.51777872*bp**(-0.80473167))
+        else if (bp .le. 20.0) then
+            per=factor*(4.73419689-0.072919483*bp)
+        else
+            per=factor*(4.73419689-0.072919483*bp)
+            print *, 'warning: pressure is too high'
+            print *, 'equation of state is only valid for < 20 GPa'
+        endif 
+        vol = 4.0*pi/3.0*br**3.0
+
+        n = vol/per
+        print *, "Adding ", n, " atoms of species ", species
+
+        !Now add n atoms randomly within the sphere 
+        do i = 1, n
+            do while(.true.)
+                do j = 1, 3
+                    call random_number(rand)
+                    p(j) = rand*(2*br) + c(j)-br
+                end do
+                if (norm2(p-c) < br) exit
+            end do
+            
+            call add_atom(0, new_type, 1, p)
+        end do
+
+    end subroutine bubble 
+
+    subroutine parse_bubble(arg_pos) 
+        integer, intent(inout) :: arg_pos
+
+        integer :: i, arglen
+        character(len=100) :: tmptxt
+
+        !First read in the bubble centroid
+        do i = 1, 3
+            arg_pos = arg_pos + 1 
+            call get_command_argument(arg_pos, tmptxt, arglen)
+            call parse_pos(i, tmptxt, c(i))
+        end do
+
+        !Now the bubble radius
+        arg_pos = arg_pos + 1
+        call get_command_argument(arg_pos, tmptxt, arglen)
+        if(arglen == 0) stop "Missing bubble radius"
+        read(tmptxt, *) br
+
+        !Now bubble pressure
+        arg_pos = arg_pos + 1
+        call get_command_argument(arg_pos, tmptxt, arglen)
+        if(arglen == 0) stop "Missing bubble pressure"
+        read(tmptxt, *) bp
+        
+        !Now bubble temperature
+        arg_pos = arg_pos + 1
+        call get_command_argument(arg_pos, tmptxt, arglen)
+        if(arglen == 0) stop "Missing bubble temperature"
+        read(tmptxt, *) bt
+
+        !Now the bubble species
+        arg_pos = arg_pos + 1
+        call get_command_argument(arg_pos, species, arglen)
+        if(arglen == 0) stop "Missing bubble species"
+
+        arg_pos = arg_pos +1 
+        return
+    end subroutine parse_bubble
+end module opt_bubble
+
diff --git a/src/opt_delete.f90 b/src/opt_delete.f90
index c3bcc9e..c111b90 100644
--- a/src/opt_delete.f90
+++ b/src/opt_delete.f90
@@ -102,7 +102,7 @@ module opt_delete
                                         for_delete(delete_num) = max(i,nei)
 
                                         !Now zero out the larger index
-                                        if(i > nei) then 
+                                        if(i < nei) then 
                                             which_cell(:,i) = 0
                                             cycle atom_loop
                                         else
diff --git a/src/parameters.f90 b/src/parameters.f90
index c589b4b..2cc6759 100644
--- a/src/parameters.f90
+++ b/src/parameters.f90
@@ -13,6 +13,7 @@ module parameters
     
     !Numeric constants
     real(kind=dp), parameter :: pi = 3.14159265358979323846_dp
+    real(kind=dp), parameter :: pvt2n =  1362.626470955822
     
     !Mode type that is being called 
     character(len=100) :: mode