Mode calc with tot_virial calculation

src/box.f90

@@ -103,4 +103,11 @@ module box
         box_bc = "ppp"
         box_bd(:) = 0
     end subroutine reset_box
+
+    pure function box_volume()
+        real(kind = dp) :: box_volume
+
+        box_volume = (box_bd(2) - box_bd(1)) * (box_bd(4) - box_bd(3))*(box_bd(6) - box_bd(5))
+        return
+    end function
 end module box

src/call_mode.f90

@@ -6,6 +6,7 @@ subroutine call_mode(arg_pos)
     use mode_convert
     use mode_merge
     use mode_metric
+    use mode_calc
     use parameters
 
     implicit none    
@@ -21,6 +22,8 @@ subroutine call_mode(arg_pos)
         call merge(arg_pos) 
     case('--metric')
         call metric(arg_pos)
+    case('--calc')
+        call calc(arg_pos)
     case default
         print *, "Mode ", trim(adjustl(mode)), " currently not accepted. Please check documentation for ", &
                  "accepted modes and rerun."

src/elements.f90

@@ -50,6 +50,9 @@ module elements
 
     !Additional module level variables we need 
     logical :: wrap_flag
+
+    !flags for data variables
+    logical :: vflag
     public 
     contains
 
@@ -136,6 +139,7 @@ module elements
         real(kind=dp), dimension(3,max_ng_node) :: adjustVar
 
         adjustVar(:,:) = 0.0_dp
+        vflag = .false.
 
         select case(trim(ele_type))
 
@@ -823,6 +827,7 @@ do i = 1, atom_num
         energy_atom(ia) = eng
         force_atom(:,ia) = force(:) 
         virial_atom(:,ia) = virial(:) 
+        vflag = .true.
         return
     end subroutine add_atom_data
 
@@ -832,6 +837,7 @@ do i = 1, atom_num
         real(kind=dp), intent(in) :: eng(max_basisnum, max_ng_node), &
                                      force(3,max_basisnum, max_ng_node), &
                                      virial(6,max_basisnum,max_ng_node)
+        vflag = .true.
         energy_node(:,:,ie)   = eng
         force_node(:,:,:,ie)  = force
         virial_node(:,:,:,ie) = virial

src/functions.f90

@@ -391,4 +391,12 @@ END FUNCTION StrDnCase
         end if
     end function permutation
 
+    pure function evtogp(virial)
+        real(kind=dp), dimension(6), intent(in) :: virial
+        real(kind=dp), dimension(6) :: evtogp
+
+        evtogp = virial * 1e21_dp * 1.602176565e-19_dp
+
+    end function
+
 end module functions

src/mode_calc.f90

@@ -0,0 +1,95 @@
+module mode_calc
+    !This mode is used to calculate various quantities based on input information 
+    use parameters 
+    use io
+    use subroutines
+    use elements
+    use box
+
+    character(len=100) :: calc_opt
+    real(kind=dp), allocatable :: calculated(:)
+    public
+    contains
+    subroutine calc(arg_pos)
+        !Main calling subroutine for mode_create
+        integer, intent(out) :: arg_pos
+
+        print *, '------------------------Mode Calc----------------------------'
+
+        !First parse command
+        call parse(arg_pos)
+        
+        print *, "Calculating ", trim(adjustl(calc_opt)), " for ", ele_num, " elements and ", atom_num, " atoms."
+        !Now call the correct calc function based on calc_opt
+        select case(trim(adjustl(calc_opt)))
+        case('tot_virial')
+            allocate(calculated(6))
+            call calc_tot_virial
+        case default
+            print *, trim(adjustl(calc_opt)), " is not accepted as a calc option in mode_calc"
+            stop 3
+        end select
+    end subroutine calc
+
+    subroutine parse(arg_pos)
+        !This parses the mode calc options
+        integer, intent(out) :: arg_pos
+
+        character(len = 100) :: infile
+        integer:: arglen
+        real(kind=dp) :: temp_box_bd(6)
+
+        call get_command_argument(2, infile, arglen)
+        if (arglen == 0 ) stop "Missing calc option in mode calc"
+        call get_in_file(infile)
+        call read_in(1, (/0.0_dp, 0.0_dp, 0.0_dp /), temp_box_bd)
+        call grow_box(temp_box_bd) 
+
+        call get_command_argument(3, calc_opt, arglen)
+        if (arglen == 0 ) stop "Missing calc option in mode calc"
+        
+        arg_pos = 4
+    end subroutine parse
+
+    subroutine calc_tot_virial
+        !Calculate the the total box pressure in GPa
+        
+        integer :: i, j, ibasis, inod
+        real(kind=dp) :: avg_virial(6)
+
+        !First check to make sure that the virial was set for the atoms/elements
+        if(.not.vflag) then 
+            print *, "Virial data has not been sent/may not be available with your current input file "
+            stop 3
+        end if
+
+        !Sum the atom virials
+        calculated = 0
+        do i = 1, atom_num
+            do j = 1, 6
+                calculated(j) = calculated(j) + virial_atom(j, i)
+            end do
+        end do
+
+        !Sum the nodal virials
+        do i = 1, ele_num
+            avg_virial(:) = 0
+            do inod = 1, ng_node(lat_ele(i))
+                do ibasis = 1, basisnum(lat_ele(i))
+                    do j = 1,6
+                        avg_virial(j) = avg_virial(j) + virial_node(j,ibasis,inod,i)/(basisnum(lat_ele(i))*ng_node(lat_ele(i)))
+                    end do
+                end do
+            end do
+
+            !Now add the total virial from the element 
+            calculated = calculated + avg_virial*(esize**3.0_dp)
+        end do
+
+        !Now calculate the total box virial and convert to GPa
+        calculated = evtogp(calculated)/box_volume()
+
+        print *, "Total virial is calculated as : (v11, v22, v33, v32, v31, v21)"
+        print *, calculated
+    end subroutine 
+end module mode_calc