support providing pbc as min and max, add test

src/fparc.f90

@@ -50,7 +50,7 @@ double precision function fparc(icart,firstjcart)
       !
       ! Otherwise, compute distance and evaluate function for this pair
       !
-      vdiff = delta_vector(xcart(icart,:), xcart(jcart,:), pbc_box)
+      vdiff = delta_vector(xcart(icart,:), xcart(jcart,:), pbc_length)
       datom = ( vdiff(1) )**2 + ( vdiff(2) )**2 + ( vdiff(3) )**2
       tol = (radius(icart)+radius(jcart))**2
       if ( datom < tol ) then

src/getinp.f90

@@ -163,8 +163,20 @@ subroutine getinp()
          read(keyword(i,3),*,iostat=ioerr) pbc_box(2)
          read(keyword(i,4),*,iostat=ioerr) pbc_box(3)
          if ( ioerr /= 0 ) exit
+         if (maxkeywords > 4 .and. trim(adjustl(keyword(i,5))) /= "none") then 
+            read(keyword(i,5),*,iostat=ioerr) pbc_box(4)
+            read(keyword(i,6),*,iostat=ioerr) pbc_box(5)
+            read(keyword(i,7),*,iostat=ioerr) pbc_box(6)
+            if ( ioerr /= 0 ) exit
+         else 
+            pbc_box(4:6) = pbc_box(1:3)
+            pbc_box(1:3) = 0.0
+         end if
+         if ( ioerr /= 0 ) exit
          using_pbc = .true.
-         write(*,*) ' Periodic boundary condition activated box: ', pbc_box(1), pbc_box(2), pbc_box(3)
+         pbc_length(1:3) = pbc_box(4:6) - pbc_box(1:3)
+         write(*,"(a, 3f8.2)") '  Periodic boundary condition activated: ', pbc_length(1), pbc_length(2), pbc_length(3)
+         write(*,"(a, 6f8.2)") '  PBC Reference box: ', pbc_box(2), pbc_box(2), pbc_box(3), pbc_box(4), pbc_box(5), pbc_box(6)
       else if( keyword(i,1) /= 'tolerance' .and. &
          keyword(i,1) /= 'short_tol_dist' .and. &
          keyword(i,1) /= 'short_tol_scale' .and. &
@@ -737,15 +749,15 @@ subroutine getinp()
       irest = irest + 1
       irestline(irest) = -1
       ityperest(irest) = 3
-      restpars(irest,1) = -dism/2
-      restpars(irest,2) = -dism/2
-      restpars(irest,3) = -dism/2
-      restpars(irest,4) = pbc_box(1) + dism/2
-      restpars(irest,5) = pbc_box(2) + dism/2
-      restpars(irest,6) = pbc_box(3) + dism/2
-      write(*,*) " For periodic boundary condition."
-      write(*,*) " We automatically add a constraint for all atoms."
-      write(*,*) " inside box -tol/2 -tol/2 -tol/2 box1+tol/2 box2+tol/2 box3+tol/2"
+      restpars(irest,1) = pbc_box(1) - dism/2
+      restpars(irest,2) = pbc_box(2) - dism/2
+      restpars(irest,3) = pbc_box(3) - dism/2
+      restpars(irest,4) = pbc_box(4) + dism/2
+      restpars(irest,5) = pbc_box(5) + dism/2
+      restpars(irest,6) = pbc_box(6) + dism/2
+      write(*,*) " PBC on: We automatically add a constraint for non-fixed atoms:"
+      write(*,"(a, 6f8.2)") "  -> inside box ", restpars(irest,1), restpars(irest,2), restpars(irest,3),&
+         restpars(irest,4), restpars(irest,5), restpars(irest,6)
    end if
 
    nrest = irest

src/gparc.f90

@@ -50,7 +50,7 @@ subroutine gparc(icart,firstjcart)
       ! Otherwise, compute distance and evaluate function for this pair
       !
       tol = (radius(icart)+radius(jcart))**2
-      vdiff = delta_vector(xcart(icart,:), xcart(jcart,:), pbc_box)
+      vdiff = delta_vector(xcart(icart,:), xcart(jcart,:), pbc_length)
       datom = (vdiff(1))**2 + (vdiff(2))**2 + (vdiff(3))**2
       if( datom < tol ) then
          dtemp = fscale(icart)*fscale(jcart) * 4.d0 * (datom - tol)

src/initial.f90

@@ -271,7 +271,7 @@ subroutine initial(n,x)
 
    if (using_pbc) then
       sizemin = 0.0d0
-      sizemax = pbc_box
+      sizemax = pbc_length
    end if
 
    ! Computing the size of the patches

src/pbc.f90

@@ -4,19 +4,19 @@
 
 module pbc
 
-   double precision, public :: pbc_box(3)
+   double precision, public :: pbc_box(6), pbc_length(3)
    logical, public :: using_pbc = .false.
 
    public delta_vector, idx_box
 
 contains
 
-   elemental double precision function delta_vector(v1,v2,pbc_box)
+   elemental double precision function delta_vector(v1,v2,pbc_length)
       implicit none
-      double precision, intent(in) :: v1, v2, pbc_box
+      double precision, intent(in) :: v1, v2, pbc_length
       delta_vector = v1 - v2
       if (using_pbc) then
-         delta_vector = delta_vector - pbc_box * nint(delta_vector/pbc_box)
+         delta_vector = delta_vector - pbc_length * nint(delta_vector/pbc_length)
       end if
    end function delta_vector
 

src/setibox.f90

@@ -15,9 +15,9 @@ subroutine setibox(x,y,z,sizemin,boxl,nboxes,iboxx,iboxy,iboxz)
    integer :: nboxes(3), iboxx, iboxy, iboxz
 
    if (using_pbc) then
-      x = x - floor(x / pbc_box(1)) * pbc_box(1)
-      y = y - floor(y / pbc_box(2)) * pbc_box(2)
-      z = z - floor(z / pbc_box(3)) * pbc_box(3)
+      x = x - floor(x / pbc_length(1)) * pbc_length(1)
+      y = y - floor(y / pbc_length(2)) * pbc_length(2)
+      z = z - floor(z / pbc_length(3)) * pbc_length(3)
    end if
 
    xtemp = x - sizemin(1)

testing/input_files/spherical_pbc.inp

@@ -0,0 +1,42 @@
+#
+# Double layered vesicle with water inside and outside
+#
+
+tolerance 2.0 
+pbc -47.5 -47.5 -47.5 47.5 47.5 47.5
+
+structure ./structure_files/palmitoil.pdb 
+  resnumbers 2
+  number 20
+  atoms 37
+    inside sphere 0. 0. 0. 14. 
+  end atoms
+  atoms 5
+    outside sphere 0. 0. 0. 26.
+  end atoms
+end structure 
+
+structure ./structure_files/palmitoil.pdb 
+  resnumbers 2
+  number 30
+  atoms 5
+    inside sphere 0. 0. 0. 29. 
+  end atoms
+  atoms 37
+    outside sphere 0. 0. 0. 41.
+  end atoms
+end structure 
+
+structure ./structure_files/water.pdb
+   resnumbers 2
+   number 50
+   inside sphere 0. 0. 0. 13.
+end structure
+
+structure ./structure_files/water.pdb
+   resnumbers 2
+   number 5000
+   outside sphere 0. 0. 0. 43.
+end structure
+
+output output.pdb

testing/runtests.jl

@@ -31,8 +31,16 @@ function check_mind(input_file::String)
         keyword, values... = split(line)
         keyword == "tolerance" && (tolerance = parse(Float64, values[1]))
         keyword == "output" && (output_name = values[1])
-        keyword == "pbc" && (unitcell = parse.(Float64,values[1:3]))
         keyword == "precision" && (precision = values[1])
+        if keyword == "pbc" 
+            if length(values) == 3
+                unitcell = parse.(Float64,values[1:3])
+            elseif length(values) == 6
+                unitcell = parse.(Float64,values[4:6]) - parse.(Float64,values[1:3])
+            else
+                error("pbc not properly set")
+            end
+        end
     end
     if isnothing(tolerance) || isnothing(output_name)
         error("tolerance or output not found")
@@ -53,14 +61,15 @@ function check_mind(input_file::String)
             
         """)
     end
-    printstyled("Test $input_file OK. \n"; color=:green, bold=true)
+    printstyled(" OK. \n"; color=:green, bold=true)
     return nothing
 end 
 
 println("Running tests...")
 if !isinteractive()
     packmol = joinpath(@__DIR__,"..","packmol")
     for input_test in ARGS
+        print(" Running test $input_test ...")
         log = IOBuffer()
         run(pipeline(`$packmol`; stdin=input_test, stdout=log))
         if occursin("Success!", String(take!(log)))

testing/test.sh

@@ -16,7 +16,8 @@ julia runtests.jl ./input_files/water_box.inp \
                   ./input_files/water_box_pbc.inp \
                   ./input_files/water_box_pbc2.inp \
                   ./input_files/bilayer_pbc.inp \
-                  ./input_files/solvprotein_pbc.inp
+                  ./input_files/solvprotein_pbc.inp \
+                  ./input_files/spherical_pbc.inp \
 
 # check if output files are properly generated in a failed run
 ./test_failed.sh ./input_files/water_box_failed.inp packmol.log "FORCED"