Change bubble code from pressure to ratio

src/opt_bubble.f90

@@ -1,6 +1,7 @@
 module opt_bubble
     !This module contains the bubble option which can be used to create voids with specific pressures of certain atoms
 
+    use atoms
     use parameters
     use elements 
     use box
@@ -8,7 +9,7 @@ module opt_bubble
 
     implicit none
 
-    real(kind=dp), private :: br, bt, bp, c(3)
+    real(kind=dp), private :: br, brat, c(3)
     character(len=2), private :: species
 
     public 
@@ -16,8 +17,8 @@ module opt_bubble
     subroutine bubble(arg_pos)
         integer, intent(inout) :: arg_pos
 
-        integer :: new_type, n, j, i
+        integer :: new_type, n, j, i, atom_num_pre
-        real(kind = dp)  :: p(3), rand, factor, per, vol
+        real(kind = dp)  :: p(3), rand, factor, per, vol, mass
 
         print *, '------------------------------------------------------------'
         print *, 'Option Bubble'
@@ -28,34 +29,36 @@ module opt_bubble
         !Now we use the existing group code to delete a sphere which represents the bubble
         centroid=c
         radius = br
-        type='both'
+        type='all'
         gshape='sphere'
         group_nodes = .true.
 
         call get_group
         call refine_group
         call get_group
+        atom_num_pre= atom_num
         call delete_group
 
         !Now we create a new atom type with the desired species
-        call add_atom_type(species, new_type)
+        call atommass(species, mass)
+        call add_atom_type(mass, 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
+        !print *, "Creating a bubble with pressure", bp, " at temperature ", bt, " with radius ", br
-        
+        !
-        factor = 1.0e24/6.02214e23
+        !factor = 1.0e24/6.02214e23
-        if (bp <= 10.0) then
+        !if (bp <= 10.0) then
-            per=factor*(3.29674113+4.51777872*bp**(-0.80473167))
+        !    per=factor*(3.29674113+4.51777872*bp**(-0.80473167))
-        else if (bp .le. 20.0) then
+        !else if (bp .le. 20.0) then
-            per=factor*(4.73419689-0.072919483*bp)
+        !    per=factor*(4.73419689-0.072919483*bp)
-        else
+        !else
-            per=factor*(4.73419689-0.072919483*bp)
+        !    per=factor*(4.73419689-0.072919483*bp)
-            print *, 'warning: pressure is too high'
+        !    print *, 'warning: pressure is too high'
-            print *, 'equation of state is only valid for < 20 GPa'
+        !    print *, 'equation of state is only valid for < 20 GPa'
-        endif 
+        !endif 
-        vol = 4.0*pi/3.0*br**3.0
+        !vol = 4.0*pi/3.0*br**3.0
-
+
-        n = vol/per
+        n = brat*(atom_num_pre-atom_num)
         print *, "Adding ", n, " atoms of species ", species
 
         !Now add n atoms randomly within the sphere 
@@ -79,34 +82,33 @@ module opt_bubble
         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)
+            print *,  trim(adjustl(tmptxt))
             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)
+        print *,  trim(adjustl(tmptxt))
         if(arglen == 0) stop "Missing bubble radius"
         read(tmptxt, *) br
 
-        !Now bubble pressure
+        !Now bubble ratio
         arg_pos = arg_pos + 1
         call get_command_argument(arg_pos, tmptxt, arglen)
-        if(arglen == 0) stop "Missing bubble pressure"
+        print *,  trim(adjustl(tmptxt))
-        read(tmptxt, *) bp
+        if(arglen == 0) stop "Missing bubble ratio"
+        read(tmptxt, *) brat
         
-        !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)
+        print *,  species
         if(arglen == 0) stop "Missing bubble species"
 
         arg_pos = arg_pos +1