Fix dangling pointers when running multiple FF calculations in a row.

Author: NikoOinonen

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

@@ -1,6 +1,7 @@
 #!/usr/bin/python
 
 import sys
+import weakref
 
 import numpy as np
 
@@ -191,6 +192,9 @@ def perform_relaxation(
         PPdisp = None
     if verbose > 0:
         print("<<<END: perform_relaxation()")
+
+    core.deleteFF_Fpointer()
+
     return fzs, PPpos, PPdisp, lvecScan
 
 
@@ -204,11 +208,14 @@ 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)
+    weakref.finalize(FF, lambda: core.deleteFF_Fpointer())  # Set array pointer to NULL when garbage collector runs.
     if Vpot:
         V = np.zeros((gridN[2], gridN[1], gridN[0]))
         core.setFF_Epointer(V)
+        weakref.finalize(FF, lambda: core.deleteFF_Epointer())
     else:
         V = None
     return FF, V
@@ -455,3 +462,5 @@ def subtractCoreDensities(
         print("sum(RHO), Nelec: ", rho.sum(), rho.sum() * dV)  # check sum
     if bSaveDebugDens:
         io.saveXSF("rho_subCoreChg.xsf", rho, lvec, head=head)
+
+    core.deleteFF_Epointer()

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

@@ -87,6 +87,24 @@ def setFF_Epointer(gridE):
     lib.setFF_Epointer(gridE)
 
 
+# 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):
     n_ = None
     if gridF is not 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()