Fix dangling pointers when running multiple FF calculations in a row (…

…#322)

The force field arrays were left as dangling pointers when running multiple computational steps in a row, which manifested in #316. This fixes the issue by setting some cleanup code to run on garbage collection. It's not very pretty since you also have to manually run the garbage collector to make sure that the pointers are deleted in time, but it works.

examples/pyridineDensOverlap/example_pyridine_density_overlap.py

@@ -32,7 +32,7 @@ def example_pyridine_density_overlap():
 
     generate_conv_rho("-s CHGCAR.xsf -t tip/density_CO.xsf -B 1.0 -E".split())
     generate_elff("-i LOCPOT.xsf --tip_dens tip/density_CO.xsf --Rcore 0.7 -E --doDensity".split())
-    generate_dftd3("-i LOCPOT.xsf --df_name PBE --energy".split())
+    generate_dftd3("-i LOCPOT.xsf --df_name PBE".split())
 
     relaxed_scan("-k 0.25 -q 1.0 --noLJ --Apauli 18.0 --bDebugFFtot".split())  # Note the --noLJ for loading separate Pauli and vdW instead of LJ force field
     plot_results("-k 0.25 -q 1.0 -a 2.0 --df".split())

ppafm/HighLevel.py

@@ -191,6 +191,9 @@ def perform_relaxation(
         PPdisp = None
     if verbose > 0:
         print("<<<END: perform_relaxation()")
+
+    core.deleteFF_Fpointer()
+
     return fzs, PPpos, PPdisp, lvecScan
 
 
@@ -204,7 +207,8 @@ def prepareArrays(FF, Vpot, parameters):
         gridN = parameters.gridN
         FF = np.zeros((gridN[2], gridN[1], gridN[0], 3))
     else:
-        parameters.gridN = np.shape(FF)
+        gridN = np.shape(FF)
+        parameters.gridN = gridN
     core.setFF_Fpointer(FF)
     if Vpot:
         V = np.zeros((gridN[2], gridN[1], gridN[0]))

ppafm/cli/conv_rho.py

@@ -1,4 +1,6 @@
 #!/usr/bin/python
+import gc
+
 import numpy as np
 
 from .. import common, fieldFFT, io
@@ -69,6 +71,10 @@ def main(argv=None):
     force_field = io.packVecGrid(f_x * args.Apauli, f_y * args.Apauli, f_z * args.Apauli)
     io.save_vec_field("FF" + namestr, force_field, lvec_sample, data_format=args.output_format, head=head_sample)
 
+    # Make sure that the energy and force field pointers are deleted so that they don't interfere if any other force fields are computed after this.
+    del energy, force_field
+    gc.collect()
+
 
 if __name__ == "__main__":
     main()

ppafm/cli/generateDFTD3.py

@@ -1,5 +1,6 @@
 #!/usr/bin/env python3
 
+import gc
 import sys
 
 from ppafm.defaults import d3
@@ -50,6 +51,9 @@ def main(argv=None):
 
     computeDFTD3(args.input, df_params=df_params, geometry_format=args.input_format, save_format=args.output_format, compute_energy=args.energy, parameters=parameters)
 
+    # Make sure that the energy and force field pointers are deleted so that they don't interfere if any other force fields are computed after this.
+    gc.collect()
+
 
 if __name__ == "__main__":
     main()

ppafm/cli/generateElFF.py

@@ -1,4 +1,5 @@
 #!/usr/bin/python
+import gc
 import sys
 from pathlib import Path
 
@@ -161,6 +162,10 @@ def main(argv=None):
     if args.energy:
         io.save_scal_field("Eel", e_electrostatic, lvec_samp, data_format=args.output_format, head=head_samp, atomic_info=(atoms_samp[:4], lvec_samp))
 
+    # Make sure that the energy and force field pointers are deleted so that they don't interfere if any other force fields are computed after this.
+    del e_electrostatic, ff_electrostatic
+    gc.collect()
+
 
 if __name__ == "__main__":
     main()

ppafm/cli/generateElFF_point_charges.py

@@ -1,5 +1,7 @@
 #!/usr/bin/python
 
+import gc
+
 from .. import common
 from ..HighLevel import computeELFF_pointCharge
 
@@ -13,6 +15,9 @@ def main(argv=None):
 
     computeELFF_pointCharge(args.input, geometry_format=args.input_format, tip=args.tip, save_format=args.output_format, computeVpot=args.energy, parameters=parameters)
 
+    # Make sure that the energy and force field pointers are deleted so that they don't interfere if any other force fields are computed after this.
+    gc.collect()
+
 
 if __name__ == "__main__":
     main()

ppafm/cli/generateLJFF.py

@@ -1,4 +1,5 @@
 #!/usr/bin/python
+import gc
 from pathlib import Path
 
 from .. import common
@@ -22,6 +23,9 @@ def main(argv=None):
         parameters=parameters,
     )
 
+    # Make sure that the energy and force field pointers are deleted so that they don't interfere if any other force fields are computed after this.
+    gc.collect()
+
 
 if __name__ == "__main__":
     main()

ppafm/core.py

@@ -1,5 +1,6 @@
 #!/usr/bin/python
 
+import weakref
 from ctypes import POINTER, Structure, c_double, c_int
 
 import numpy as np
@@ -76,6 +77,7 @@ def setFF_shape(n_, cell, parameters=None):
 
 def setFF_Fpointer(gridF):
     lib.setFF_Fpointer(gridF)
+    weakref.finalize(gridF, deleteFF_Fpointer)  # Set array pointer to NULL when garbage collector runs.
 
 
 # void setFF_pointer( double * gridF, double * gridE  )
@@ -85,6 +87,25 @@ def setFF_Fpointer(gridF):
 
 def setFF_Epointer(gridE):
     lib.setFF_Epointer(gridE)
+    weakref.finalize(gridE, deleteFF_Epointer)  # Set array pointer to NULL when garbage collector runs.
+
+
+# void deleteFF_Fpointer()
+lib.deleteFF_Fpointer.argtypes = []
+lib.deleteFF_Fpointer.restype = None
+
+
+def deleteFF_Fpointer():
+    lib.deleteFF_Fpointer()
+
+
+# void deleteFF_Epointer()
+lib.deleteFF_Epointer.argtypes = []
+lib.deleteFF_Epointer.restype = None
+
+
+def deleteFF_Epointer():
+    lib.deleteFF_Epointer()
 
 
 def setFF(cell=None, gridF=None, gridE=None, parameters=None):

ppafm/cpp/ProbeParticle.cpp

@@ -344,6 +344,16 @@ DLLEXPORT void setFF_Epointer( double * gridE_ ){
     gridE = gridE_;
 }
 
+// set force field array pointer to NULL
+DLLEXPORT void deleteFF_Fpointer(){
+    gridF = NULL;
+}
+
+// set energy array pointer to NULL
+DLLEXPORT void deleteFF_Epointer(){
+    gridE = NULL;
+}
+
 // set forcefield grid dimension "n"
 DLLEXPORT void setGridN( int * n ){
     //gridShape.n.set( *(Vec3i*)n );

tests/_test_vdw.py

@@ -4,6 +4,8 @@
 Compare the C++ and OpenCL implementations of the vdW calculation and check that they are consistent.
 """
 
+import gc
+
 import numpy as np
 import pyopencl as cl
 
@@ -102,3 +104,6 @@ def test_dftd3():
     assert np.allclose(coeffs_ocl, coeffs_cpp)
     assert np.allclose(FF_ocl[..., :3], FF_cpp, rtol=1e-4, atol=1e-6)
     assert np.allclose(FF_ocl[..., 3], E_cpp, rtol=1e-4, atol=1e-6)
+
+    del E_cpp, FF_cpp
+    gc.collect()