Fix pointers

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()