Merge pull request #30 from aselimov/ft--update-boundary-command

Ft  update boundary command

src/Makefile

@@ -1,6 +1,6 @@
 FC=ifort
-FFLAGS=-mcmodel=large -g -O0 -stand f08 -fpe0 -traceback -check bounds,uninit -warn all -implicitnone -no-wrap-margin -heap-arrays
+#FFLAGS=-mcmodel=large -g -O0 -stand f08 -fpe0 -traceback -check bounds,uninit -warn all -implicitnone -no-wrap-margin -heap-arrays
-#FFLAGS=-mcmodel=large -Ofast -no-wrap-margin -heap-arrays
+FFLAGS=-mcmodel=large -Ofast -no-wrap-margin -heap-arrays
 MODES=mode_create.o mode_merge.o mode_convert.o
 OPTIONS=opt_disl.o opt_group.o opt_orient.o opt_delete.o 
 OBJECTS=main.o elements.o io.o subroutines.o functions.o atoms.o call_mode.o box.o  $(MODES) $(OPTIONS) call_option.o

src/box.f90

@@ -6,6 +6,7 @@ module box
 
     real(kind=dp) :: box_bd(6) !Global box boundaries
     character(len=3) :: box_bc !Box boundary conditions (periodic or shrinkwrapped)
+    logical :: bound_called 
     !The subbox variables contain values for each subbox, being the boxes read in through some
     !command. Currently only mode_merge will require sub_boxes, for mode_create it will always 
     !allocate to only 1 sub_box
@@ -27,6 +28,7 @@ module box
         !Initialize some box functions
         box_bd(:) = 0.0_dp
         box_bc = 'ppp'
+        bound_called=.false.
     end subroutine box_init
 
     subroutine alloc_sub_box(n)

src/call_option.f90

@@ -28,6 +28,7 @@ subroutine call_option(option, arg_pos)
         arg_pos=arg_pos+1
         call get_command_argument(arg_pos, box_bc)
         arg_pos=arg_pos+1
+        bound_called = .true.
     case('-delete')
         call run_delete(arg_pos)
     case default

src/elements.f90

@@ -468,8 +468,8 @@ module elements
 
         do i = 1, atom_num
             do j = 1, 3
-                if (r_atom(j,i) > max_bd(j)) max_bd(j) = r_atom(j,i) + lim_zero
+                if (r_atom(j,i) > max_bd(j)) max_bd(j) = r_atom(j,i) + tol 
-                if (r_atom(j,i) < min_bd(j)) min_bd(j) = r_atom(j,i) - lim_zero
+                if (r_atom(j,i) < min_bd(j)) min_bd(j) = r_atom(j,i) - tol 
             end do 
         end do 
 
@@ -477,17 +477,18 @@ module elements
             do inod = 1, ng_node(lat_ele(i))
                 do ibasis = 1, basisnum(lat_ele(i))
                     do j = 1, 3
-                        if (r_node(j,ibasis,inod,i) > max_bd(j)) max_bd(j) = r_node(j,ibasis,inod,i) + lim_zero
+                        if (r_node(j,ibasis,inod,i) > max_bd(j)) max_bd(j) = r_node(j,ibasis,inod,i) + tol
-                        if (r_node(j,ibasis,inod,i) < min_bd(j)) min_bd(j) = r_node(j,ibasis,inod,i) -lim_zero
+                        if (r_node(j,ibasis,inod,i) < min_bd(j)) min_bd(j) = r_node(j,ibasis,inod,i) - tol
-
                     end do 
                 end do 
             end do 
         end do 
 
         do j = 1, 3
-            box_bd(2*j) = max_bd(j)
+            if(box_bc(j:j) == 's') then 
-            box_bd(2*j-1) = min_bd(j)
+                box_bd(2*j) = max_bd(j)
+                box_bd(2*j-1) = min_bd(j)
+            end if
         end do 
     end subroutine
 

src/main.f90

@@ -101,6 +101,9 @@ program main
     !If wrap flag was passed then call the wrap atoms command
     if(wrap_flag) call wrap_atoms
 
+    !If we called the boundary command then we adjust the box bounds
+    if(bound_called) call def_new_box
+
     !Check to make sure a file was passed to be written out and then write out 
             ! Before building do a check on the file
     if (outfilenum == 0) then

src/opt_orient.f90

@@ -22,6 +22,7 @@ module opt_orient
         integer :: i, ibasis, inod
         logical :: isortho, isrighthanded
         real(kind=dp) :: inv_sub_box_ori(3,3,sub_box_num)
+        character(len=3) :: old_box_bc
 
         !First parse the orient command
         call parse_orient(arg_pos)
@@ -31,7 +32,7 @@ module opt_orient
 
         !Find all inverse orientation matrices for all sub_boxes
         do i = 1, sub_box_num
-            call matrix_inverse(sub_box_ori, 3, inv_sub_box_ori)
+            call matrix_inverse(sub_box_ori(:,:,i), 3, inv_sub_box_ori(:,:,i))
         end do
 
         !Now transform all atoms
@@ -62,8 +63,11 @@ module opt_orient
         !Save original box boundaries
         orig_box_bd = box_bd
 
-        !Now find new box boundaries
+        !Now find new box boundaries, have to temporarily define the box as shrink wrapped for def new box to work
+        old_box_bc = box_Bc
+        box_bc = 'sss'
         call def_new_box
+        box_bc = old_box_bc
     end subroutine orient
 
     subroutine parse_orient(arg_pos)

src/parameters.f90

@@ -6,7 +6,8 @@ module parameters
     integer, parameter :: dp= selected_real_kind(15,307)
     !Parameters for floating point tolerance
     real(kind=dp), parameter :: lim_zero = epsilon(1.0_dp), &
-                                lim_large = huge(1.0_dp)
+                                lim_large = huge(1.0_dp), &
+                                tol = 10.0_dp**(-6.0_dp)
     logical, save :: lmpcac 
     
     !Numeric constants

src/subroutines.f90

@@ -230,11 +230,13 @@ module subroutines
         integer :: j 
         real(kind=dp) ::box_len
         do j = 1, 3
-            box_len = box_bd(2*j) - box_bd(2*j-1)
+            if(box_bc(j:j) == 'p') then 
-            if (r(j) > box_bd(2*j)) then 
+                box_len = box_bd(2*j) - box_bd(2*j-1)
-                r(j) = r(j) - box_len
+                if (r(j) > box_bd(2*j)) then 
-            else if (r(j) < box_bd(2*j-1)) then 
+                    r(j) = r(j) - box_len
-                r(j) = r(j) + box_len
+                else if (r(j) < box_bd(2*j-1)) then 
+                    r(j) = r(j) + box_len
+                end if
             end if
         end do
     end subroutine