diff --git a/src/box.f90 b/src/box.f90 index 83d81a8..c179136 100644 --- a/src/box.f90 +++ b/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 diff --git a/src/call_mode.f90 b/src/call_mode.f90 index 8d73e25..b933751 100644 --- a/src/call_mode.f90 +++ b/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." diff --git a/src/elements.f90 b/src/elements.f90 index 62b56f0..5d3b3c6 100644 --- a/src/elements.f90 +++ b/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 diff --git a/src/functions.f90 b/src/functions.f90 index 0a8bc49..ebd9a50 100644 --- a/src/functions.f90 +++ b/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 diff --git a/src/mode_calc.f90 b/src/mode_calc.f90 new file mode 100644 index 0000000..497c90b --- /dev/null +++ b/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