From 326fb2adf75e69e9e1d785216149fadd9efc81d0 Mon Sep 17 00:00:00 2001 From: Martin Ondracek Date: Thu, 22 May 2025 17:57:01 +0200 Subject: [PATCH 1/9] Rearange array indices and memory layout in CLI so that vectors fields are stored as F[ix,iy,iz,i] and scalar fields as V[ix,iy,iz], both in C layout --- ppafm/GUIWidgets.py | 12 ++-- ppafm/GridUtils.py | 6 +- ppafm/HighLevel.py | 39 ++++++------- ppafm/PPPlot.py | 107 +++++++++++++++++++++++++++--------- ppafm/cli/generateElFF.py | 8 +-- ppafm/cli/plot_results.py | 18 +++--- ppafm/common.py | 16 +++--- ppafm/cpp/Grid.h | 98 ++++++++++++++++++++------------- ppafm/cpp/GridUtils.cpp | 14 ++--- ppafm/cpp/ProbeParticle.cpp | 76 ++++++++++++------------- ppafm/cpp/integerOps.h | 2 +- ppafm/fieldFFT.py | 20 +++---- ppafm/io.py | 20 +++---- ppafm/ocl/AFMulator.py | 3 +- ppafm/ocl/field.py | 4 +- 15 files changed, 254 insertions(+), 189 deletions(-) diff --git a/ppafm/GUIWidgets.py b/ppafm/GUIWidgets.py index ae6745df..b8d7cb4b 100644 --- a/ppafm/GUIWidgets.py +++ b/ppafm/GUIWidgets.py @@ -173,10 +173,10 @@ def plotSlice(self, F_stack, z_slice, title=None, margins=None, grid_selector=0, if margins: self.axes.add_patch( matplotlib.patches.Rectangle( - (margins[0], margins[1]), F.shape[1] - margins[2] - margins[0], F.shape[0] - margins[3] - margins[1], linewidth=2, edgecolor="r", facecolor="none" + (margins[0], margins[1]), F.shape[0] - margins[2] - margins[0], F.shape[1] - margins[3] - margins[1], linewidth=2, edgecolor="r", facecolor="none" ) ) - textRes = "output size: " + str(F.shape[1] - margins[2] - margins[0]) + "x" + str(F.shape[0] - margins[3] - margins[1]) + textRes = "output size: " + str(F.shape[0] - margins[2] - margins[0]) + "x" + str(F.shape[1] - margins[3] - margins[1]) if slice_length: textRes += " length [A] =" + f"{slice_length[0]:03.4f}, {slice_length[1]:03.4f}" self.axes.set_xlabel(textRes) @@ -202,8 +202,8 @@ def plotSlice2(self, F_stack, title=None, margins=None, grid_selector=0, slice_l F = F * (1 - alpha) + big_len_image * alpha self.img = self.axes.imshow(F, origin="lower", cmap="viridis", interpolation="bicubic") - j_min, i_min = np.unravel_index(F.argmin(), F.shape) - j_max, i_max = np.unravel_index(F.argmax(), F.shape) + i_min, j_min = np.unravel_index(F.argmin(), F.shape) + i_max, j_max = np.unravel_index(F.argmax(), F.shape) if margins: self.axes.add_patch(matplotlib.patches.Rectangle((margins[0], margins[1]), margins[2], margins[3], linewidth=2, edgecolor="r", facecolor="none")) @@ -212,8 +212,8 @@ def plotSlice2(self, F_stack, title=None, margins=None, grid_selector=0, slice_l textRes += " length [A] =" + f"{slice_length[0]:03.4f}, {slice_length[1]:03.4f}" self.axes.set_xlabel(textRes) - self.axes.set_xlim(0, F.shape[1]) - self.axes.set_ylim(0, F.shape[0]) + self.axes.set_xlim(0, F.shape[0]) + self.axes.set_ylim(0, F.shape[1]) self.axes.set_title(title) if grid_selector > 0: self.axes.grid(True, linestyle="dotted", color="blue") diff --git a/ppafm/GridUtils.py b/ppafm/GridUtils.py index 845eab59..875492fd 100644 --- a/ppafm/GridUtils.py +++ b/ppafm/GridUtils.py @@ -156,6 +156,6 @@ def sphericalHist(data, center, dr, n): def readNumsUpTo(filename, dimensions, noline): - N_arry = np.zeros((dimensions[0] * dimensions[1] * dimensions[2]), dtype=np.double) - lib.ReadNumsUpTo_C(filename.encode(), N_arry, dimensions, noline) - return N_arry + N_array = np.zeros((dimensions[0] * dimensions[1] * dimensions[2]), dtype=np.double) + lib.ReadNumsUpTo_C(filename.encode(), N_array, dimensions, noline) + return N_array diff --git a/ppafm/HighLevel.py b/ppafm/HighLevel.py index f150bd8f..92fd294a 100755 --- a/ppafm/HighLevel.py +++ b/ppafm/HighLevel.py @@ -68,9 +68,9 @@ def relaxedScan3D(xTips, yTips, zTips, trj=None, bF3d=False): rTips = np.zeros((ntips, 3)) rs = np.zeros((ntips, 3)) fs = np.zeros((ntips, 3)) - nx = len(zTips) + nx = len(xTips) ny = len(yTips) - nz = len(xTips) + nz = len(zTips) if bF3d: fzs = np.zeros((nx, ny, nz, 3)) else: @@ -90,14 +90,14 @@ def relaxedScan3D(xTips, yTips, zTips, trj=None, bF3d=False): rTips[:, 2] = trj[:, 2] core.relaxTipStroke(rTips, rs, fs) if bF3d: - fzs[:, iy, ix, 0] = (fs[:, 0].copy())[::-1] - fzs[:, iy, ix, 1] = (fs[:, 1].copy())[::-1] - fzs[:, iy, ix, 2] = (fs[:, 2].copy())[::-1] + fzs[ix, iy, :, 0] = fs[::-1, 0] + fzs[ix, iy, :, 1] = fs[::-1, 1] + fzs[ix, iy, :, 2] = fs[::-1, 2] else: - fzs[:, iy, ix] = (fs[:, 2].copy())[::-1] - PPpos[:, iy, ix, 0] = rs[::-1, 0] - PPpos[:, iy, ix, 1] = rs[::-1, 1] - PPpos[:, iy, ix, 2] = rs[::-1, 2] + fzs[ix, iy, :] = fs[::-1, 2] + PPpos[ix, iy, :, 0] = rs[::-1, 0] + PPpos[ix, iy, :, 1] = rs[::-1, 1] + PPpos[ix, iy, :, 2] = rs[::-1, 2] if verbose > 0: print("<<>BEGIN: relaxedScan3D_omp()") if verbose > 0: print(" zTips : ", zTips) - nz = len(zTips) - ny = len(yTips) nx = len(xTips) + ny = len(yTips) + nz = len(zTips) rTips = np.zeros((nx, ny, nz, 3)) rs = np.zeros((nx, ny, nz, 3)) fs = np.zeros((nx, ny, nz, 3)) @@ -118,14 +118,11 @@ def relaxedScan3D_omp(xTips, yTips, zTips, trj=None, bF3d=False, tip_spline=None rTips[:, :, :, 1] = yTips[None, :, None] rTips[:, :, :, 2] = zTips[::-1][None, None, :] core.relaxTipStrokes_omp(rTips, rs, fs, tip_spline=tip_spline) - # rs[:,:,:,:] = rs[:,:,::-1,:].transpose(2,1,0,3).copy() - # fs[:,:,:,:] = fs[:,:,::-1,:] - # fzs = fs[:,:,::-1,2].transpose(2,1,0).copy() - rs = rs[:, :, ::-1, :].transpose(2, 1, 0, 3).copy() + rs = rs[:, :, ::-1, :].copy() if bF3d: - fzs = fs[:, :, ::-1, :].transpose(2, 1, 0, 3).copy() + fzs = fs[:, :, ::-1, :].copy() else: - fzs = fs[:, :, ::-1, 2].transpose(2, 1, 0).copy() + fzs = fs[:, :, ::-1, 2].copy() if verbose > 0: print("<<0 zpos[i] %= lvec[3, 2] - ff_kpfm_t0sv[i, :, :] = ff_kpfm_t0sv[i, :, :] / ((v_ref_s) * (zpos[i] + 0.1)) - ff_kpfm_tvs0[i, :, :] = ff_kpfm_tvs0[i, :, :] / ((v_ref_t) * (zpos[i] + 0.1)) + ff_kpfm_t0sv[:, :, i] = ff_kpfm_t0sv[:, :, i] / ((v_ref_s) * (zpos[i] + 0.1)) + ff_kpfm_tvs0[:, :, i] = ff_kpfm_tvs0[:, :, i] / ((v_ref_t) * (zpos[i] + 0.1)) print(">>> Saving electrostatic forcefield ... ") io.save_vec_field("FFkpfm_t0sV", ff_kpfm_t0sv, lvec_samp, data_format=args.output_format, head=head_samp) diff --git a/ppafm/cli/plot_results.py b/ppafm/cli/plot_results.py index fc973775..6a97f53b 100644 --- a/ppafm/cli/plot_results.py +++ b/ppafm/cli/plot_results.py @@ -141,7 +141,7 @@ def main(argv=None): dirname + "/xy" + atoms_str + cbar_str, pp_positions[:, :, :, 0], pp_positions[:, :, :, 1], - slices=list(range(0, len(pp_positions))), + slices=list(range(0, pp_positions.shape[2])), BG=pp_positions[:, :, :, 2], extent=extent, atoms=atoms, @@ -164,7 +164,7 @@ def main(argv=None): PPPlot.plotImages( dirname + "/IETS" + atoms_str + cbar_str, iets, - slices=list(range(0, len(iets))), + slices=list(range(0, iets.shape[2])), zs=tip_positions_z, extent=extent, atoms=atoms, @@ -177,7 +177,7 @@ def main(argv=None): PPPlot.plotImages( dirname + "/Evib" + atoms_str + cbar_str, e_vib[:, :, :, 0], - slices=list(range(0, len(iets))), + slices=list(range(0, iets.shape[2])), zs=tip_positions_z, extent=extent, atoms=atoms, @@ -190,7 +190,7 @@ def main(argv=None): PPPlot.plotImages( dirname + "/Kvib" + atoms_str + cbar_str, 16.0217662 * eigenvalue_k[:, :, :, 0], - slices=list(range(0, len(iets))), + slices=(range(0, iets.shape[2])), zs=tip_positions_z, extent=extent, atoms=atoms, @@ -206,7 +206,7 @@ def main(argv=None): PPPlot.plotImages( dirname + "/OutI" + atoms_str + cbar_str, current, - slices=list(range(0, len(current))), + slices=(range(0, current.shape[2])), zs=tip_positions_z, extent=extent, atoms=atoms, @@ -344,7 +344,7 @@ def main(argv=None): PPPlot.plotImages( dir_name_amplitude + "/df" + atoms_str + cbar_str, dfs, - slices=list(range(0, len(dfs))), + slices=list(range(0, dfs.shape[2])), zs=tip_positions_z + parameters.Amplitude / 2.0, extent=extent, cmap=parameters.colorscale, @@ -405,7 +405,7 @@ def main(argv=None): PPPlot.plotImages( dir_name_lcpd + "/LCPD" + atoms_str + cbar_str, lcpd, - slices=list(range(0, len(lcpd))), + slices=list(range(0, lcpd.shape[2])), zs=tip_positions_z + parameters.Amplitude / 2.0, extent=extent, cmap=parameters.colorscale_kpfm, @@ -422,7 +422,7 @@ def main(argv=None): PPPlot.plotImages( dir_name_lcpd + "/_Asym-LCPD" + atoms_str + cbar_str, lcpd, - slices=list(range(0, len(lcpd))), + slices=list(range(0, lcpd.shape[2])), zs=tip_positions_z + parameters.Amplitude / 2.0, extent=extent, cmap=parameters.colorscale_kpfm, @@ -459,7 +459,7 @@ def main(argv=None): PPPlot.plotImages( dirname + "/Fz" + atoms_str + cbar_str, fzs, - slices=list(range(0, len(fzs))), + slices=listlist(range(0, fzs.shape[2])), zs=tip_positions_z, # + parameters.Amplitude / 2.0, # no oscillation to the force extent=extent, cmap=parameters.colorscale, diff --git a/ppafm/common.py b/ppafm/common.py index ed5b7ec0..86d1ec3c 100644 --- a/ppafm/common.py +++ b/ppafm/common.py @@ -433,7 +433,7 @@ def Fz2df(F, dz, k0, f0, amplitude=1.0, units=16.0217656): Internal force units are eV/A the 16.021... converts stifness eV/A**2 to N/m """ W = get_df_weight(amplitude, dz=dz) - dFconv = np.apply_along_axis(lambda m: np.convolve(m, W, mode="valid"), axis=0, arr=F) + dFconv = np.apply_along_axis(lambda m: np.convolve(m, W, mode="valid"), axis=2, arr=F) return dFconv * units * f0 / k0 @@ -446,9 +446,9 @@ def Fz2df_tilt(F, d, k0, f0, amplitude=1.0, units=16.0217656): """ dr = np.sqrt(d[0] ** 2 + d[1] ** 2 + d[2] ** 2) W = get_df_weight(amplitude, dz=dr) - dFconv_x = np.apply_along_axis(lambda m: np.convolve(m, W, mode="valid"), axis=0, arr=F[:, :, :, 0]) - dFconv_y = np.apply_along_axis(lambda m: np.convolve(m, W, mode="valid"), axis=0, arr=F[:, :, :, 1]) - dFconv_z = np.apply_along_axis(lambda m: np.convolve(m, W, mode="valid"), axis=0, arr=F[:, :, :, 2]) + dFconv_x = np.apply_along_axis(lambda m: np.convolve(m, W, mode="valid"), axis=2, arr=F[:, :, :, 0]) + dFconv_y = np.apply_along_axis(lambda m: np.convolve(m, W, mode="valid"), axis=2, arr=F[:, :, :, 1]) + dFconv_z = np.apply_along_axis(lambda m: np.convolve(m, W, mode="valid"), axis=2, arr=F[:, :, :, 2]) return (dFconv_x * d[0] + dFconv_y * d[1] + dFconv_z * d[2]) * units * f0 / k0 @@ -1000,11 +1000,11 @@ def genFFSampling(lvec, pixPerAngstrome=10): def getPos(lvec, nDim=None, pixPerAngstrome=10): if nDim is None: nDim = genFFSampling(lvec, pixPerAngstrome=pixPerAngstrome) - dCell = np.array((lvec[1, :] / nDim[2], lvec[2, :] / nDim[1], lvec[3, :] / nDim[0])) + dCell = np.array((lvec[1, :] / nDim[0], lvec[2, :] / nDim[1], lvec[3, :] / nDim[2])) ABC = np.mgrid[0 : nDim[0], 0 : nDim[1], 0 : nDim[2]] - X = lvec[0, 0] + ABC[2] * dCell[0, 0] + ABC[1] * dCell[1, 0] + ABC[0] * dCell[2, 0] - Y = lvec[0, 1] + ABC[2] * dCell[0, 1] + ABC[1] * dCell[1, 1] + ABC[0] * dCell[2, 1] - Z = lvec[0, 2] + ABC[2] * dCell[0, 2] + ABC[1] * dCell[1, 2] + ABC[0] * dCell[2, 2] + X = lvec[0, 0] + ABC[0] * dCell[0, 0] + ABC[1] * dCell[1, 0] + ABC[2] * dCell[2, 0] + Y = lvec[0, 1] + ABC[0] * dCell[0, 1] + ABC[1] * dCell[1, 1] + ABC[2] * dCell[2, 1] + Z = lvec[0, 2] + ABC[0] * dCell[0, 2] + ABC[1] * dCell[1, 2] + ABC[2] * dCell[2, 2] return X, Y, Z diff --git a/ppafm/cpp/Grid.h b/ppafm/cpp/Grid.h index 80b3e9c4..78101494 100644 --- a/ppafm/cpp/Grid.h +++ b/ppafm/cpp/Grid.h @@ -15,7 +15,7 @@ #define fast_floor_offset 1000 #define fast_floor( x ) ( ( (int)( x + fast_floor_offset ) ) - fast_floor_offset ) -#define i3D( ix, iy, iz ) ( iz*nxy + iy*nx + ix ) +#define i3D( ix, iy, iz ) ( (ix*ny + iy)*nz + iz ) // ================= CONSTANTS @@ -73,11 +73,21 @@ inline double interpolate3DWrap( double * grid, const Vec3i& n, const Vec3d& r ) int xoff = n.x<<3; int imx = r.x +xoff; double tx = r.x - imx +xoff; double mx = 1 - tx; int itx = (imx+1)%n.x; imx=imx%n.x; int yoff = n.y<<3; int imy = r.y +yoff; double ty = r.y - imy +yoff; double my = 1 - ty; int ity = (imy+1)%n.y; imy=imy%n.y; int zoff = n.z<<3; int imz = r.z +zoff; double tz = r.z - imz +zoff; double mz = 1 - tz; int itz = (imz+1)%n.z; imz=imz%n.z; - int nxy = n.x * n.y; int nx = n.x; - //double out = grid[ i3D( imx, imy, imz ) ]; - - //double out = mz * my * ( ( mx * grid[ i3D( imx, imy, imz ) ] ) + ( tx * grid[ i3D( itx, imy, imz ) ] ); - //ty * ( mx * grid[ i3D( imx, ity, imz ) ] ) + ( tx * grid[ i3D( itx, ity, imz ) ] ) ); + int nx = n.x; int ny = n.y; int nz = n.z; + + /* + int imx, imy, imz, itx, ity, itz; + doulble mx, my, mz, tx, ty, tz; + if(r.x >= 0) {imx = r.x; tx = r.x - imx; mx = 1 - tx; imx = imx % nx;} + else {itx = r.x; mx = itx - r.x; tx = 1 - mx; imx = itx % nx + n.x-1;} + itx = (imx+1) % nx; + if(r.y >= 0) {imy = r.y; ty = r.y - imy; my = 1 - ty; imy = imy % ny;} + else {ity = r.y; my = ity - r.y; ty = 1 - my; imy = ity % ny + n.y-1;} + ity = (imy+1) % ny; + if(r.z >= 0) {imz = r.z; tz = r.z - imz; mz = 1 - tz; imz = imz % nz;} + else {itz = r.z; mz = itz - r.z; tz = 1 - mz; imz = itz % nz + n.z-1;} + itz = (imz+1) % nz; + */ double out = mz * ( my * ( ( mx * grid[ i3D( imx, imy, imz ) ] ) + ( tx * grid[ i3D( itx, imy, imz ) ] ) ) + @@ -96,7 +106,22 @@ inline Vec3d interpolate3DvecWrap( Vec3d * grid, const Vec3i& n, const Vec3d& r int xoff = n.x<<3; int imx = r.x +xoff; double tx = r.x - imx +xoff; double mx = 1 - tx; int itx = (imx+1)%n.x; imx=imx%n.x; int yoff = n.y<<3; int imy = r.y +yoff; double ty = r.y - imy +yoff; double my = 1 - ty; int ity = (imy+1)%n.y; imy=imy%n.y; int zoff = n.z<<3; int imz = r.z +zoff; double tz = r.z - imz +zoff; double mz = 1 - tz; int itz = (imz+1)%n.z; imz=imz%n.z; - int nxy = n.x * n.y; int nx = n.x; + + int nx = n.x; int ny = n.y; int nz = n.z; + /* + int imx, imy, imz, itx, ity, itz; + doulble mx, my, mz, tx, ty, tz; + if(r.x >= 0) {imx = r.x; tx = r.x - imx; mx = 1 - tx; imx = imx % nx;} + else {itx = r.x; mx = itx - r.x; tx = 1 - mx; imx = itx % nx + n.x-1;} + itx = (imx+1) % nx; + if(r.y >= 0) {imy = r.y; ty = r.y - imy; my = 1 - ty; imy = imy % ny;} + else {ity = r.y; my = ity - r.y; ty = 1 - my; imy = ity % ny + n.y-1;} + ity = (imy+1) % ny; + if(r.z >= 0) {imz = r.z; tz = r.z - imz; mz = 1 - tz; imz = imz % nz;} + else {itz = r.z; mz = itz - r.z; tz = 1 - mz; imz = itz % nz + n.z-1;} + itz = (imz+1) % nz; + */ + double mymx = my*mx; double mytx = my*tx; double tymx = ty*mx; double tytx = ty*tx; Vec3d out; out.set_mul( grid[ i3D( imx, imy, imz ) ], mz*mymx ); out.add_mul( grid[ i3D( itx, imy, imz ) ], mz*mytx ); @@ -110,48 +135,45 @@ inline Vec3d interpolate3DvecWrap( Vec3d * grid, const Vec3i& n, const Vec3d& r // iterate over field template< void FUNC( int ibuff, const Vec3d& pos_, void * args ) > -void interateGrid3D( const Vec3d& pos0, const Vec3i& n, const Mat3d& dCell, void * args ){ +void iterateGrid3D( const Vec3d& pos0, const Vec3i& n, const Mat3d& dCell, void * args ){ int nx = n.x; int ny = n.y; int nz = n.z; - //int nx = n.z; int ny = n.y; int nz = n.x; - int nxy = ny * nx; - printf( "interateGrid3D nx,y,z (%i,%i,%i) nxy %i\n", nx,ny,nz, nxy ); + printf( "iterateGrid3D nx,y,z (%i,%i,%i)\n", nx, ny, nz); Vec3d pos; pos.set( pos0 ); - //printf(" interateGrid3D : args %i \n", args ); - for ( int ic=0; ic -void interateGrid3D_omp( const Vec3d& pos0, const Vec3i& n, const Mat3d& dCell, void * args ){ +void iterateGrid3D_omp( const Vec3d& pos0, const Vec3i& n, const Mat3d& dCell, void * args ){ int ntot = n.x*n.y*n.z; - int ncpu = omp_get_num_threads(); printf( "interateGrid3D_omp nx,y,z (%i,%i,%i) nxy %i ncpu %i \n", n.x,n.y,n.z, ntot, ncpu ); + int ncpu = omp_get_num_threads(); printf( "iterateGrid3D_omp nx,y,z (%i,%i,%i) ntot %i ncpu %i \n", n.x,n.y,n.z, ntot, ncpu ); int ndone=0; #pragma omp parallel for collapse(3) shared(pos0,n,dCell,args,ndone) - for ( int ic=0; ic( r0, gridShape.n, gridShape.dCell, NULL ); + iterateGrid3D( r0, gridShape.n, gridShape.dCell, NULL ); } // --------- find center of mass @@ -165,7 +161,7 @@ extern "C" { DLLEXPORT double cog( double * data_, double* center ){ data = data_; Histogram::Htot += 0; Histogram::center.set(0.0); Vec3d r0; r0.set(0.0,0.0,0.0); - interateGrid3D( r0, gridShape.n, gridShape.dCell, NULL ); + iterateGrid3D( r0, gridShape.n, gridShape.dCell, NULL ); Histogram::center.mul( 1/Histogram::Htot ); ((Vec3d*)center)->set(Histogram::center); return Histogram::Htot; @@ -181,7 +177,7 @@ extern "C" { DLLEXPORT void setGridN( int * n ){ //gridShape.n.set( *(Vec3i*)n ); - gridShape.n.set( n[2], n[1], n[0] ); + gridShape.n.set( n[0], n[1], n[2] ); printf( " nxyz %i %i %i \n", gridShape.n.x, gridShape.n.y, gridShape.n.z ); } diff --git a/ppafm/cpp/ProbeParticle.cpp b/ppafm/cpp/ProbeParticle.cpp index b75c968b..d284a94f 100644 --- a/ppafm/cpp/ProbeParticle.cpp +++ b/ppafm/cpp/ProbeParticle.cpp @@ -357,7 +357,7 @@ DLLEXPORT void deleteFF_Epointer(){ // set forcefield grid dimension "n" DLLEXPORT void setGridN( int * n ){ //gridShape.n.set( *(Vec3i*)n ); - gridShape.n.set( n[2], n[1], n[0] ); + gridShape.n.set( n[0], n[1], n[2] ); printf( " nxyz %i %i %i \n", gridShape.n.x, gridShape.n.y, gridShape.n.z ); } @@ -508,31 +508,31 @@ DLLEXPORT void computeD3Coeffs( DLLEXPORT void getLennardJonesFF( int natoms_, double * Ratoms_, double * cLJs ){ natoms=natoms_; Ratoms=(Vec3d*)Ratoms_; nCoefPerAtom = 2; Vec3d r0; r0.set(0.0,0.0,0.0); - //interateGrid3D < evalCell < addAtom_LJ > >( r0, gridShape.n, gridShape.dCell, cLJs ); - interateGrid3D_omp < evalCell < addAtom_LJ > >( r0, gridShape.n, gridShape.dCell, cLJs ); + //iterateGrid3D < evalCell < addAtom_LJ > >( r0, gridShape.n, gridShape.dCell, cLJs ); + iterateGrid3D_omp < evalCell < addAtom_LJ > >( r0, gridShape.n, gridShape.dCell, cLJs ); } DLLEXPORT void getVdWFF( int natoms_, double * Ratoms_, double * cLJs ){ natoms=natoms_; Ratoms=(Vec3d*)Ratoms_; nCoefPerAtom = 2; Vec3d r0; r0.set(0.0,0.0,0.0); - //interateGrid3D < evalCell < addAtom_VdW > >( r0, gridShape.n, gridShape.dCell, cLJs ); - interateGrid3D_omp < evalCell < addAtom_VdW > >( r0, gridShape.n, gridShape.dCell, cLJs ); + //iterateGrid3D < evalCell < addAtom_VdW > >( r0, gridShape.n, gridShape.dCell, cLJs ); + iterateGrid3D_omp < evalCell < addAtom_VdW > >( r0, gridShape.n, gridShape.dCell, cLJs ); } DLLEXPORT void getDFTD3FF(int natoms_, double * Ratoms_, double *d3_coeffs){ natoms = natoms_; Ratoms = (Vec3d*)Ratoms_; nCoefPerAtom = 4; Vec3d r0; r0.set(0.0,0.0,0.0); - //interateGrid3D < evalCell < addAtom_DFTD3 > >( r0, gridShape.n, gridShape.dCell, d3_coeffs ); - interateGrid3D_omp < evalCell < addAtom_DFTD3 > >( r0, gridShape.n, gridShape.dCell, d3_coeffs ); + //iterateGrid3D < evalCell < addAtom_DFTD3 > >( r0, gridShape.n, gridShape.dCell, d3_coeffs ); + iterateGrid3D_omp < evalCell < addAtom_DFTD3 > >( r0, gridShape.n, gridShape.dCell, d3_coeffs ); } DLLEXPORT void getMorseFF( int natoms_, double * Ratoms_, double * REs, double alpha ){ natoms=natoms_; Ratoms=(Vec3d*)Ratoms_; nCoefPerAtom = 2; Morse_alpha = alpha; Vec3d r0; r0.set(0.0,0.0,0.0); - //interateGrid3D < evalCell < addAtom_Morse > >( r0, gridShape.n, gridShape.dCell, REs ); - interateGrid3D_omp < evalCell < addAtom_Morse > >( r0, gridShape.n, gridShape.dCell, REs ); + //iterateGrid3D < evalCell < addAtom_Morse > >( r0, gridShape.n, gridShape.dCell, REs ); + iterateGrid3D_omp < evalCell < addAtom_Morse > >( r0, gridShape.n, gridShape.dCell, REs ); } // sample Coulomb Force-field on 3D mesh over provided set of atoms with positions Rs_[i] with constant kQQs = - k_coulomb * Q_ProbeParticle * Q[i] @@ -542,14 +542,14 @@ DLLEXPORT void getCoulombFF( int natoms_, double * Ratoms_, double * kQQs, int k Vec3d r0; r0.set(0.0,0.0,0.0); //printf(" kind %i \n", kind ); switch(kind){ - //case 0: interateGrid3D < evalCell < foo > >( r0, gridShape.n, gridShape.dCell, kQQs_ ); - //case 0: interateGrid3D < evalCell < addAtom_Coulomb_s > >( r0, gridShape.n, gridShape.dCell, kQQs ); break; - //case 1: interateGrid3D < evalCell < addAtom_Coulomb_pz > >( r0, gridShape.n, gridShape.dCell, kQQs ); break; - //case 2: interateGrid3D < evalCell < addAtom_Coulomb_dz2 > >( r0, gridShape.n, gridShape.dCell, kQQs ); break; - - case 0: interateGrid3D_omp < evalCell < addAtom_Coulomb_s > >( r0, gridShape.n, gridShape.dCell, kQQs ); break; - case 1: interateGrid3D_omp < evalCell < addAtom_Coulomb_pz > >( r0, gridShape.n, gridShape.dCell, kQQs ); break; - case 2: interateGrid3D_omp < evalCell < addAtom_Coulomb_dz2 > >( r0, gridShape.n, gridShape.dCell, kQQs ); break; + //case 0: iterateGrid3D < evalCell < foo > >( r0, gridShape.n, gridShape.dCell, kQQs_ ); + //case 0: iterateGrid3D < evalCell < addAtom_Coulomb_s > >( r0, gridShape.n, gridShape.dCell, kQQs ); break; + //case 1: iterateGrid3D < evalCell < addAtom_Coulomb_pz > >( r0, gridShape.n, gridShape.dCell, kQQs ); break; + //case 2: iterateGrid3D < evalCell < addAtom_Coulomb_dz2 > >( r0, gridShape.n, gridShape.dCell, kQQs ); break; + + case 0: iterateGrid3D_omp < evalCell < addAtom_Coulomb_s > >( r0, gridShape.n, gridShape.dCell, kQQs ); break; + case 1: iterateGrid3D_omp < evalCell < addAtom_Coulomb_pz > >( r0, gridShape.n, gridShape.dCell, kQQs ); break; + case 2: iterateGrid3D_omp < evalCell < addAtom_Coulomb_dz2 > >( r0, gridShape.n, gridShape.dCell, kQQs ); break; } } @@ -581,20 +581,20 @@ DLLEXPORT void getVdWFF_RE( int natoms_, double * Ratoms_, double * REs, int kin //if(ADamp>0){ ADamp = ADamp_; } switch(kind){ //case 0: if(ADamp_>0){ ADamp_Const = ADamp_; }; E=addAtom_VdW ( dR, fout, coefs ); break; - //case 1: if(ADamp_>0){ ADamp_R2 = ADamp_; } interateGrid3D>( r0, gridShape.n, gridShape.dCell, REs ); break; - //case 2: if(ADamp_>0){ ADamp_R4 = ADamp_; } interateGrid3D>( r0, gridShape.n, gridShape.dCell, REs ); break; - //case 3: if(ADamp_>0){ ADamp_invR4 = ADamp_; } interateGrid3D>( r0, gridShape.n, gridShape.dCell, REs ); break; - //case 4: if(ADamp_>0){ ADamp_invR8 = ADamp_; } interateGrid3D>( r0, gridShape.n, gridShape.dCell, REs ); break; - - case 1: if(ADamp_>0){ ADamp_R2 = ADamp_; } interateGrid3D_omp>( r0, gridShape.n, gridShape.dCell, REs ); break; - case 2: if(ADamp_>0){ ADamp_R4 = ADamp_; } interateGrid3D_omp>( r0, gridShape.n, gridShape.dCell, REs ); break; - case 3: if(ADamp_>0){ ADamp_invR4 = ADamp_; } interateGrid3D_omp>( r0, gridShape.n, gridShape.dCell, REs ); break; - case 4: if(ADamp_>0){ ADamp_invR8 = ADamp_; } interateGrid3D_omp>( r0, gridShape.n, gridShape.dCell, REs ); break; - - // case 0: interateGrid3D >>>( r0, gridShape.n, gridShape.dCell, REs ); break; - // case 1: interateGrid3D >>>( r0, gridShape.n, gridShape.dCell, REs ); break; - // case 2: interateGrid3D>>( r0, gridShape.n, gridShape.dCell, REs ); break; - // case 2: interateGrid3D>>( r0, gridShape.n, gridShape.dCell, REs ); break; + //case 1: if(ADamp_>0){ ADamp_R2 = ADamp_; } iterateGrid3D>( r0, gridShape.n, gridShape.dCell, REs ); break; + //case 2: if(ADamp_>0){ ADamp_R4 = ADamp_; } iterateGrid3D>( r0, gridShape.n, gridShape.dCell, REs ); break; + //case 3: if(ADamp_>0){ ADamp_invR4 = ADamp_; } iterateGrid3D>( r0, gridShape.n, gridShape.dCell, REs ); break; + //case 4: if(ADamp_>0){ ADamp_invR8 = ADamp_; } iterateGrid3D>( r0, gridShape.n, gridShape.dCell, REs ); break; + + case 1: if(ADamp_>0){ ADamp_R2 = ADamp_; } iterateGrid3D_omp>( r0, gridShape.n, gridShape.dCell, REs ); break; + case 2: if(ADamp_>0){ ADamp_R4 = ADamp_; } iterateGrid3D_omp>( r0, gridShape.n, gridShape.dCell, REs ); break; + case 3: if(ADamp_>0){ ADamp_invR4 = ADamp_; } iterateGrid3D_omp>( r0, gridShape.n, gridShape.dCell, REs ); break; + case 4: if(ADamp_>0){ ADamp_invR8 = ADamp_; } iterateGrid3D_omp>( r0, gridShape.n, gridShape.dCell, REs ); break; + + // case 0: iterateGrid3D >>>( r0, gridShape.n, gridShape.dCell, REs ); break; + // case 1: iterateGrid3D >>>( r0, gridShape.n, gridShape.dCell, REs ); break; + // case 2: iterateGrid3D>>( r0, gridShape.n, gridShape.dCell, REs ); break; + // case 2: iterateGrid3D>>( r0, gridShape.n, gridShape.dCell, REs ); break; } } @@ -602,8 +602,8 @@ DLLEXPORT void getGaussDensity( int natoms_, double * Ratoms_, double * cRAs ){ natoms=natoms_; Ratoms=(Vec3d*)Ratoms_; nCoefPerAtom = 2; Vec3d r0; r0.set(0.0,0.0,0.0); Vec3d* gridF_=gridF; gridF=0; - //interateGrid3D < evalCell < addAtom_Gauss > >( r0, gridShape.n, gridShape.dCell, cRAs ); - interateGrid3D_omp < evalCell < addAtom_Gauss > >( r0, gridShape.n, gridShape.dCell, cRAs ); + //iterateGrid3D < evalCell < addAtom_Gauss > >( r0, gridShape.n, gridShape.dCell, cRAs ); + iterateGrid3D_omp < evalCell < addAtom_Gauss > >( r0, gridShape.n, gridShape.dCell, cRAs ); gridF=gridF_; } @@ -611,8 +611,8 @@ DLLEXPORT void getSlaterDensity( int natoms_, double * Ratoms_, double * cRAs ){ natoms=natoms_; Ratoms=(Vec3d*)Ratoms_; nCoefPerAtom = 2; Vec3d r0; r0.set(0.0,0.0,0.0); Vec3d* gridF_=gridF; gridF=0; - //interateGrid3D < evalCell < addAtom_Slater > >( r0, gridShape.n, gridShape.dCell, cRAs ); - interateGrid3D_omp < evalCell < addAtom_Slater > >( r0, gridShape.n, gridShape.dCell, cRAs ); + //iterateGrid3D < evalCell < addAtom_Slater > >( r0, gridShape.n, gridShape.dCell, cRAs ); + iterateGrid3D_omp < evalCell < addAtom_Slater > >( r0, gridShape.n, gridShape.dCell, cRAs ); gridF=gridF_; } @@ -620,8 +620,8 @@ DLLEXPORT void getDensityR4spline( int natoms_, double * Ratoms_, double * cRAs natoms=natoms_; Ratoms=(Vec3d*)Ratoms_; nCoefPerAtom = 2; Vec3d r0; r0.set(0.0,0.0,0.0); Vec3d* gridF_=gridF; gridF=0; - //interateGrid3D < evalCell < addAtom_splineR4 > >( r0, gridShape.n, gridShape.dCell, cRAs ); - interateGrid3D_omp < evalCell < addAtom_splineR4 > >( r0, gridShape.n, gridShape.dCell, cRAs ); + //iterateGrid3D < evalCell < addAtom_splineR4 > >( r0, gridShape.n, gridShape.dCell, cRAs ); + iterateGrid3D_omp < evalCell < addAtom_splineR4 > >( r0, gridShape.n, gridShape.dCell, cRAs ); gridF=gridF_; } @@ -684,7 +684,7 @@ DLLEXPORT int relaxTipStrokes_omp( int nx, int ny, int probeStart, int relaxAlg, #pragma omp parallel for collapse(2) shared( nx, ny, probeStart, relaxAlg, nstep, rTips_, rs_, fs_, ndone ) for (int ix=0; ix -// ========= index poerations and modulo-math =========== +// ========= index operations and modulo-math =========== inline int wrap_index_fast( int i, int n){ if(i<0){ return n+i; }else if (i>=n){ return i-n; }; diff --git a/ppafm/fieldFFT.py b/ppafm/fieldFFT.py index b49cb2d2..8709a268 100644 --- a/ppafm/fieldFFT.py +++ b/ppafm/fieldFFT.py @@ -31,20 +31,20 @@ def getSize(inp_axis, dims, sampleSize): def getMGrid(dims, dd): "returns coordinate arrays X, Y, Z" (dx, dy, dz) = dd - nDim = [dims[2], dims[1], dims[0]] + nDim = [dims[0], dims[1], dims[2]] XYZ = np.mgrid[0 : nDim[0], 0 : nDim[1], 0 : nDim[2]].astype(float) # fmt: off - xshift = nDim[2]//2; xshift_ = xshift; + xshift = nDim[0]//2; xshift_ = xshift; yshift = nDim[1]//2; yshift_ = yshift; - zshift = nDim[0]//2; zshift_ = zshift; + zshift = nDim[2]//2; zshift_ = zshift; - if( nDim[2]%2 != 0 ): xshift_ += 1.0 + if( nDim[0]%2 != 0 ): xshift_ += 1.0 if( nDim[1]%2 != 0 ): yshift_ += 1.0 - if( nDim[0]%2 != 0 ): zshift_ += 1.0 + if( nDim[2]%2 != 0 ): zshift_ += 1.0 - X = dx*np.roll( XYZ[2] - xshift_, xshift, axis=2) + X = dx*np.roll( XYZ[0] - xshift_, xshift, axis=2) Y = dy*np.roll( XYZ[1] - yshift_, yshift, axis=1) - Z = dz*np.roll( XYZ[0] - zshift_, zshift, axis=0) + Z = dz*np.roll( XYZ[2] - zshift_, zshift, axis=0) # fmt: on return X, Y, Z @@ -142,7 +142,7 @@ def addCoreDensity(rho, val, x, y, z, sigma=0.1): def addCoreDensities(atoms, valElDict, rho, lvec, sigma=0.1): sampleSize = getSampleDimensions(lvec) nDim = rho.shape - dims = (nDim[2], nDim[1], nDim[0]) + dims = (nDim[0], nDim[1], nDim[2]) xsize, dx = getSize("x", dims, sampleSize) ysize, dy = getSize("y", dims, sampleSize) zsize, dz = getSize("z", dims, sampleSize) @@ -261,7 +261,7 @@ def potential2forces(V, lvec, nDim, sigma=0.7, rho=None, multipole=None, tilt=0. if verbose > 0: print("--- Preprocessing ---") sampleSize = getSampleDimensions(lvec) - dims = (nDim[2], nDim[1], nDim[0]) + dims = (nDim[0], nDim[1], nDim[2]) xsize, dx = getSize("x", dims, sampleSize) ysize, dy = getSize("y", dims, sampleSize) zsize, dz = getSize("z", dims, sampleSize) @@ -287,7 +287,7 @@ def potential2forces(V, lvec, nDim, sigma=0.7, rho=None, multipole=None, tilt=0. def potential2forces_mem(V, lvec, nDim, sigma=0.7, rho=None, multipole=None, doForce=True, doPot=False, deleteV=True, tilt=0.0): print("--- Preprocessing ---") sampleSize = getSampleDimensions(lvec) - dims = (nDim[2], nDim[1], nDim[0]) + dims = (nDim[0], nDim[1], nDim[2]) xsize, dx = getSize("x", dims, sampleSize) ysize, dy = getSize("y", dims, sampleSize) diff --git a/ppafm/io.py b/ppafm/io.py index 79763cb1..66c9660a 100644 --- a/ppafm/io.py +++ b/ppafm/io.py @@ -594,7 +594,7 @@ def saveXSF(fname, data, lvec=None, dd=None, head=XSF_HEAD_DEFAULT, verbose=1): for line in head: fileout.write(line) nDim = np.shape(data) - _writeArr(fileout, (nDim[2] + 1, nDim[1] + 1, nDim[0] + 1)) + _writeArr(fileout, (nDim[0] + 1, nDim[1] + 1, nDim[2] + 1)) _writeArr2D(fileout, lvec) data2 = np.zeros(np.array(nDim) + 1) # These crazy 3 lines are here since the first and the last cube in XSF in every direction is the same @@ -602,17 +602,16 @@ def saveXSF(fname, data, lvec=None, dd=None, head=XSF_HEAD_DEFAULT, verbose=1): data2[-1, :, :] = data2[0, :, :] data2[:, -1, :] = data2[:, 0, :] data2[:, :, -1] = data2[:, :, 0] - for r in data2.flat: + for r in data2.flatten(order="F"): fileout.write("%10.5e\n" % r) fileout.write(" END_DATAGRID_3D\n") fileout.write("END_BLOCK_DATAGRID_3D\n") -def loadXSF(fname, xyz_order=False, verbose=True): +def loadXSF(fname, verbose=True): filein = open(fname) startline, head = _readUpTo(filein, "BEGIN_DATAGRID_3D") # startline - number of the line with DATAGRID_3D_. Dinensions are located in the next line nDim = [int(iii) for iii in filein.readline().split()] # reading 1 line with dimensions - nDim.reverse() nDim = np.array(nDim) lvec = _readmat(filein, 4) # reading 4 lines where 1st line is origin of datagrid and 3 next lines are the cell vectors filein.close() @@ -623,11 +622,9 @@ def loadXSF(fname, xyz_order=False, verbose=True): F = readNumsUpTo(fname, nDim.astype(np.int32).copy(), startline + 5) if verbose: print("io | Done") - FF = np.reshape(F, nDim)[:-1, :-1, :-1] - if xyz_order: - FF = FF.transpose((2, 1, 0)) + FF = np.reshape(F, nDim, order="F")[:-1, :-1, :-1] # FF is not C_CONTIGUOUS without copy - FF = FF.copy() + FF = FF.copy(order="C") return FF, lvec, nDim - 1, head @@ -651,7 +648,7 @@ def getFromHead_PRIMCOORD(head): # =================== Cube -def loadCUBE(fname, xyz_order=False, verbose=True): +def loadCUBE(fname, verbose=True): filein = open(fname) # First two lines of the header are comments filein.readline() @@ -680,11 +677,8 @@ def loadCUBE(fname, xyz_order=False, verbose=True): if verbose: print("io | nDim: ", nDim) - FF = np.reshape(F, nDim) - if not xyz_order: - FF = FF.transpose((2, 1, 0)).copy() # Transposition of the array to have the same order of data as in XSF file + FF = np.reshape(F, nDim, order="C") - nDim = [nDim[2], nDim[1], nDim[0]] # Setting up the corresponding dimensions. head = [] head.append("BEGIN_BLOCK_DATAGRID_3D \n") head.append("g98_3D_unknown \n") diff --git a/ppafm/ocl/AFMulator.py b/ppafm/ocl/AFMulator.py index 7a08c07b..02497ee6 100644 --- a/ppafm/ocl/AFMulator.py +++ b/ppafm/ocl/AFMulator.py @@ -659,6 +659,7 @@ def plot_images(self, X, outdir="afm_images", prefix="df"): zs=zTips, extent=extent, cmap=self.colorscale, + reversed_ind=True, ) @@ -844,4 +845,4 @@ def quick_afm( scan_window[1][1], ] zs = np.linspace(scan_window[0][2], scan_window[1][2], scan_dim[2] + 1)[:num_heights] - plotImages(os.path.join(out_dir, "df"), X, range(len(X)), extent=extent, zs=zs) + plotImages(os.path.join(out_dir, "df"), X, range(len(X)), extent=extent, zs=zs, reversed_ind=True) diff --git a/ppafm/ocl/field.py b/ppafm/ocl/field.py index 54e3d258..a9e244e4 100644 --- a/ppafm/ocl/field.py +++ b/ppafm/ocl/field.py @@ -227,10 +227,10 @@ def from_file(cls, file_path, scale=1.0): file_path = str(file_path) if file_path.endswith(".cube"): - data, lvec, _, _ = io.loadCUBE(file_path, xyz_order=True, verbose=False) + data, lvec, _, _ = io.loadCUBE(file_path, verbose=False) Zs, x, y, z, _ = io.loadAtomsCUBE(file_path) elif file_path.endswith(".xsf"): - data, lvec, _, _ = io.loadXSF(file_path, xyz_order=True, verbose=False) + data, lvec, _, _ = io.loadXSF(file_path, verbose=False) try: (Zs, x, y, z, _), _, _ = io.loadXSFGeom(file_path) except ValueError: From 88a2591534341b9c4e2e8f0b93a017001a7054c5 Mon Sep 17 00:00:00 2001 From: Martin Ondracek Date: Thu, 22 May 2025 18:18:22 +0200 Subject: [PATCH 2/9] Change the way indices are calculated in interpolate3DWrap and interpolate3DvecWrap, so as to prevent SegFault with the pyridineDensOverlap example. --- ppafm/cpp/Grid.h | 114 +++++++++++++++++++++++------------------------ 1 file changed, 55 insertions(+), 59 deletions(-) diff --git a/ppafm/cpp/Grid.h b/ppafm/cpp/Grid.h index 78101494..3cf2f1b2 100644 --- a/ppafm/cpp/Grid.h +++ b/ppafm/cpp/Grid.h @@ -68,69 +68,65 @@ class GridShape { // interpolation of vector force-field Vec3d[ix,iy,iz] in periodic boundary condition inline double interpolate3DWrap( double * grid, const Vec3i& n, const Vec3d& r ){ - //#pragma omp simd - //{ - int xoff = n.x<<3; int imx = r.x +xoff; double tx = r.x - imx +xoff; double mx = 1 - tx; int itx = (imx+1)%n.x; imx=imx%n.x; - int yoff = n.y<<3; int imy = r.y +yoff; double ty = r.y - imy +yoff; double my = 1 - ty; int ity = (imy+1)%n.y; imy=imy%n.y; - int zoff = n.z<<3; int imz = r.z +zoff; double tz = r.z - imz +zoff; double mz = 1 - tz; int itz = (imz+1)%n.z; imz=imz%n.z; - int nx = n.x; int ny = n.y; int nz = n.z; - - /* - int imx, imy, imz, itx, ity, itz; - doulble mx, my, mz, tx, ty, tz; - if(r.x >= 0) {imx = r.x; tx = r.x - imx; mx = 1 - tx; imx = imx % nx;} - else {itx = r.x; mx = itx - r.x; tx = 1 - mx; imx = itx % nx + n.x-1;} - itx = (imx+1) % nx; - if(r.y >= 0) {imy = r.y; ty = r.y - imy; my = 1 - ty; imy = imy % ny;} - else {ity = r.y; my = ity - r.y; ty = 1 - my; imy = ity % ny + n.y-1;} - ity = (imy+1) % ny; - if(r.z >= 0) {imz = r.z; tz = r.z - imz; mz = 1 - tz; imz = imz % nz;} - else {itz = r.z; mz = itz - r.z; tz = 1 - mz; imz = itz % nz + n.z-1;} - itz = (imz+1) % nz; - */ - - double out = mz * ( - my * ( ( mx * grid[ i3D( imx, imy, imz ) ] ) + ( tx * grid[ i3D( itx, imy, imz ) ] ) ) + - ty * ( ( mx * grid[ i3D( imx, ity, imz ) ] ) + ( tx * grid[ i3D( itx, ity, imz ) ] ) ) ) - + tz * ( - my * ( ( mx * grid[ i3D( imx, imy, itz ) ] ) + ( tx * grid[ i3D( itx, imy, itz ) ] ) ) + - ty * ( ( mx * grid[ i3D( imx, ity, itz ) ] ) + ( tx * grid[ i3D( itx, ity, itz ) ] ) ) ); - //} - return out; + //#pragma omp simd + //{ + //int xoff = n.x<<3; int imx = r.x +xoff; double tx = r.x - imx +xoff; double mx = 1 - tx; int itx = (imx+1)%n.x; imx=imx%n.x; + //int yoff = n.y<<3; int imy = r.y +yoff; double ty = r.y - imy +yoff; double my = 1 - ty; int ity = (imy+1)%n.y; imy=imy%n.y; + //int zoff = n.z<<3; int imz = r.z +zoff; double tz = r.z - imz +zoff; double mz = 1 - tz; int itz = (imz+1)%n.z; imz=imz%n.z; + int nx = n.x; int ny = n.y; int nz = n.z; + + int imx, imy, imz, itx, ity, itz; + double mx, my, mz, tx, ty, tz; + if(r.x >= 0) {imx = r.x; tx = r.x - imx; mx = 1 - tx; imx = imx % nx;} + else {itx = r.x; mx = itx - r.x; tx = 1 - mx; imx = itx % nx + n.x-1;} + itx = (imx+1) % nx; + if(r.y >= 0) {imy = r.y; ty = r.y - imy; my = 1 - ty; imy = imy % ny;} + else {ity = r.y; my = ity - r.y; ty = 1 - my; imy = ity % ny + n.y-1;} + ity = (imy+1) % ny; + if(r.z >= 0) {imz = r.z; tz = r.z - imz; mz = 1 - tz; imz = imz % nz;} + else {itz = r.z; mz = itz - r.z; tz = 1 - mz; imz = itz % nz + n.z-1;} + itz = (imz+1) % nz; + + double out = mz * ( + my * ( ( mx * grid[ i3D( imx, imy, imz ) ] ) + ( tx * grid[ i3D( itx, imy, imz ) ] ) ) + + ty * ( ( mx * grid[ i3D( imx, ity, imz ) ] ) + ( tx * grid[ i3D( itx, ity, imz ) ] ) ) ) + + tz * ( + my * ( ( mx * grid[ i3D( imx, imy, itz ) ] ) + ( tx * grid[ i3D( itx, imy, itz ) ] ) ) + + ty * ( ( mx * grid[ i3D( imx, ity, itz ) ] ) + ( tx * grid[ i3D( itx, ity, itz ) ] ) ) ); + //} + return out; } // interpolation of vector force-field Vec3d[ix,iy,iz] in periodic boundary condition inline Vec3d interpolate3DvecWrap( Vec3d * grid, const Vec3i& n, const Vec3d& r ){ - //#pragma omp simd - //{ - int xoff = n.x<<3; int imx = r.x +xoff; double tx = r.x - imx +xoff; double mx = 1 - tx; int itx = (imx+1)%n.x; imx=imx%n.x; - int yoff = n.y<<3; int imy = r.y +yoff; double ty = r.y - imy +yoff; double my = 1 - ty; int ity = (imy+1)%n.y; imy=imy%n.y; - int zoff = n.z<<3; int imz = r.z +zoff; double tz = r.z - imz +zoff; double mz = 1 - tz; int itz = (imz+1)%n.z; imz=imz%n.z; - - int nx = n.x; int ny = n.y; int nz = n.z; - /* - int imx, imy, imz, itx, ity, itz; - doulble mx, my, mz, tx, ty, tz; - if(r.x >= 0) {imx = r.x; tx = r.x - imx; mx = 1 - tx; imx = imx % nx;} - else {itx = r.x; mx = itx - r.x; tx = 1 - mx; imx = itx % nx + n.x-1;} - itx = (imx+1) % nx; - if(r.y >= 0) {imy = r.y; ty = r.y - imy; my = 1 - ty; imy = imy % ny;} - else {ity = r.y; my = ity - r.y; ty = 1 - my; imy = ity % ny + n.y-1;} - ity = (imy+1) % ny; - if(r.z >= 0) {imz = r.z; tz = r.z - imz; mz = 1 - tz; imz = imz % nz;} - else {itz = r.z; mz = itz - r.z; tz = 1 - mz; imz = itz % nz + n.z-1;} - itz = (imz+1) % nz; - */ - - double mymx = my*mx; double mytx = my*tx; double tymx = ty*mx; double tytx = ty*tx; - Vec3d out; - out.set_mul( grid[ i3D( imx, imy, imz ) ], mz*mymx ); out.add_mul( grid[ i3D( itx, imy, imz ) ], mz*mytx ); - out.add_mul( grid[ i3D( imx, ity, imz ) ], mz*tymx ); out.add_mul( grid[ i3D( itx, ity, imz ) ], mz*tytx ); - out.add_mul( grid[ i3D( imx, ity, itz ) ], tz*tymx ); out.add_mul( grid[ i3D( itx, ity, itz ) ], tz*tytx ); - out.add_mul( grid[ i3D( imx, imy, itz ) ], tz*mymx ); out.add_mul( grid[ i3D( itx, imy, itz ) ], tz*mytx ); - //printf( "DEBUG interpolate3DvecWrap gp(%g,%g,%g) igp(%i,%i,%i)/(%i,%i,%i) %i->%g out(%g,%g,%g) \n", r.x, r.y, r.z, imx, imy, imz, n.x,n.y,n.z, i3D( imx, imy, imz ), grid[ i3D( imx, imy, imz ) ], out.x,out.y,out.z ); - //} - return out; + //#pragma omp simd + //{ + //int xoff = n.x<<3; int imx = r.x +xoff; double tx = r.x - imx +xoff; double mx = 1 - tx; int itx = (imx+1)%n.x; imx=imx%n.x; + //int yoff = n.y<<3; int imy = r.y +yoff; double ty = r.y - imy +yoff; double my = 1 - ty; int ity = (imy+1)%n.y; imy=imy%n.y; + //int zoff = n.z<<3; int imz = r.z +zoff; double tz = r.z - imz +zoff; double mz = 1 - tz; int itz = (imz+1)%n.z; imz=imz%n.z; + + int nx = n.x; int ny = n.y; int nz = n.z; + int imx, imy, imz, itx, ity, itz; + double mx, my, mz, tx, ty, tz; + if(r.x >= 0) {imx = r.x; tx = r.x - imx; mx = 1 - tx; imx = imx % nx;} + else {itx = r.x; mx = itx - r.x; tx = 1 - mx; imx = itx % nx + n.x-1;} + itx = (imx+1) % nx; + if(r.y >= 0) {imy = r.y; ty = r.y - imy; my = 1 - ty; imy = imy % ny;} + else {ity = r.y; my = ity - r.y; ty = 1 - my; imy = ity % ny + n.y-1;} + ity = (imy+1) % ny; + if(r.z >= 0) {imz = r.z; tz = r.z - imz; mz = 1 - tz; imz = imz % nz;} + else {itz = r.z; mz = itz - r.z; tz = 1 - mz; imz = itz % nz + n.z-1;} + itz = (imz+1) % nz; + + double mymx = my*mx; double mytx = my*tx; double tymx = ty*mx; double tytx = ty*tx; + Vec3d out; + out.set_mul( grid[ i3D( imx, imy, imz ) ], mz*mymx ); out.add_mul( grid[ i3D( itx, imy, imz ) ], mz*mytx ); + out.add_mul( grid[ i3D( imx, ity, imz ) ], mz*tymx ); out.add_mul( grid[ i3D( itx, ity, imz ) ], mz*tytx ); + out.add_mul( grid[ i3D( imx, ity, itz ) ], tz*tymx ); out.add_mul( grid[ i3D( itx, ity, itz ) ], tz*tytx ); + out.add_mul( grid[ i3D( imx, imy, itz ) ], tz*mymx ); out.add_mul( grid[ i3D( itx, imy, itz ) ], tz*mytx ); + //printf( "DEBUG interpolate3DvecWrap gp(%g,%g,%g) igp(%i,%i,%i)/(%i,%i,%i) %i->%g out(%g,%g,%g) \n", r.x, r.y, r.z, imx, imy, imz, n.x,n.y,n.z, i3D( imx, imy, imz ), grid[ i3D( imx, imy, imz ) ], out.x,out.y,out.z ); + //} + return out; } // iterate over field From 583af51bfeb7b4c81673f8bae01b0892a2122a7e Mon Sep 17 00:00:00 2001 From: Martin Ondracek Date: Fri, 23 May 2025 14:20:55 +0200 Subject: [PATCH 3/9] Do forgotten switching of axes for np.roll inside getMGrid to complete index order redefinition --- ppafm/fieldFFT.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/ppafm/fieldFFT.py b/ppafm/fieldFFT.py index 8709a268..782f753b 100644 --- a/ppafm/fieldFFT.py +++ b/ppafm/fieldFFT.py @@ -42,9 +42,9 @@ def getMGrid(dims, dd): if( nDim[1]%2 != 0 ): yshift_ += 1.0 if( nDim[2]%2 != 0 ): zshift_ += 1.0 - X = dx*np.roll( XYZ[0] - xshift_, xshift, axis=2) + X = dx*np.roll( XYZ[0] - xshift_, xshift, axis=0) Y = dy*np.roll( XYZ[1] - yshift_, yshift, axis=1) - Z = dz*np.roll( XYZ[2] - zshift_, zshift, axis=0) + Z = dz*np.roll( XYZ[2] - zshift_, zshift, axis=2) # fmt: on return X, Y, Z From 18d93305502abbe1f2d33d2cd486ac41fe938f9b Mon Sep 17 00:00:00 2001 From: Martin Ondracek Date: Mon, 26 May 2025 14:54:56 +0200 Subject: [PATCH 4/9] Correct index order for the capacitance term in cli/plot_results.py and memory layout in saveWSxM functions in io.py --- ppafm/cli/plot_results.py | 6 +++--- ppafm/io.py | 8 ++++---- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/ppafm/cli/plot_results.py b/ppafm/cli/plot_results.py index 6a97f53b..4e79739c 100644 --- a/ppafm/cli/plot_results.py +++ b/ppafm/cli/plot_results.py @@ -224,14 +224,14 @@ def main(argv=None): else: # Prepare to calculate KPFM/LCPD # For each pixel, find a,b,c such that df = aV^2 + bV + c - # This is done as polynomial (2nd order) regression, the above equality need not be exact + # This is done as polynomial (2nd order) regression, the above equality needs not be exact # but the least-square criterion will be used. # The coefficients a,b,c are going to be determined as linar combinations of df(V) at different biases: # a = Sum_i w_KPFM_a (V_i) # b = Sum_i w_KPFM_b (V_i) # c = Sum_i w_KPFM_c (V_i) # Now, the coefficients (weights) w_KPFM have to be determined. - # This will be done with the help of Gram-sSchmidt ortogonalization: + # This will be done with the help of Gram-Schmidt ortogonalization: # Create vectors v0, v1, v2, # v0 = [1]_(i=1..N), # v1 = [V_i]_(i=1..N,) @@ -316,7 +316,7 @@ def main(argv=None): if applied_bias: r_tip = parameters.Rtip for iz, z in enumerate(tip_positions_z): - fzs[iz, :, :] = fzs[iz, :, :] - np.pi * parameters.permit * ((r_tip * r_tip) / ((z - args.z0) * (z + r_tip))) * (voltage - args.V0) * ( + fzs[:, :, iz] = fzs[:, :, iz] - np.pi * parameters.permit * ((r_tip * r_tip) / ((z - args.z0) * (z + r_tip))) * (voltage - args.V0) * ( voltage - args.V0 ) dfs = common.Fz2df( diff --git a/ppafm/io.py b/ppafm/io.py index 66c9660a..5430c968 100644 --- a/ppafm/io.py +++ b/ppafm/io.py @@ -691,10 +691,10 @@ def loadCUBE(fname, verbose=True): def saveWSxM_2D(name_file, data, Xs, Ys): - tmp_data = data.flatten() + tmp_data = data.flatten(order="F") out_data = np.zeros((len(tmp_data), 3)) - out_data[:, 0] = Xs.flatten() - out_data[:, 1] = Ys.flatten() + out_data[:, 0] = Xs.flatten(order="F") + out_data[:, 1] = Ys.flatten(order="F") out_data[:, 2] = tmp_data # .copy() f = open(name_file, "w") print("WSxM file copyright Nanotec Electronica", file=f) @@ -715,7 +715,7 @@ def saveWSxM_3D(prefix, data, extent, slices=None): for i in slices: print("slice no: ", i) fname = prefix + "_%03d.xyz" % i - saveWSxM_2D(fname, data[i], Xs, Ys) + saveWSxM_2D(fname, data[:, :, i], Xs, Ys) # ================ Npy From 1125f5c6d01b01cc0d2fcfcc0c76bf0f5234637e Mon Sep 17 00:00:00 2001 From: Martin Ondracek Date: Tue, 6 May 2025 19:17:33 +0200 Subject: [PATCH 5/9] Repair the bug in the --rotate option of ppafm-relaxed-scan by flipping the array indices of force fields --- ppafm/cli/relaxed_scan.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/ppafm/cli/relaxed_scan.py b/ppafm/cli/relaxed_scan.py index 8f886855..bd0b13a7 100755 --- a/ppafm/cli/relaxed_scan.py +++ b/ppafm/cli/relaxed_scan.py @@ -104,34 +104,34 @@ def main(argv=None): print("Loading Pauli force field from FFpauli_{x,y,z}") ff_pauli, lvec, _, atomic_info_or_head = io.load_vec_field("FFpauli", data_format=args.output_format) - ff_pauli[0, :, :, :], ff_pauli[1, :, :, :] = rotate_ff(ff_pauli[0, :, :, :], ff_pauli[1, :, :, :], opt_dict["rotate"]) + ff_pauli[:, :, :, 0], ff_pauli[1, :, :, :, 1] = rotate_ff(ff_pauli[:, :, :, 0], ff_pauli[:, :, :, 1], opt_dict["rotate"]) print("Loading vdW force field from FFvdW_{x,y,z}") ff_vdw, lvec, _, atomic_info_or_head = io.load_vec_field("FFvdW", data_format=args.output_format) - ff_vdw[0, :, :, :], ff_vdw[1, :, :, :] = rotate_ff(ff_vdw[0, :, :, :], ff_vdw[1, :, :, :], opt_dict["rotate"]) + ff_vdw[:, :, :, 0], ff_vdw[:, :, :, 1] = rotate_ff(ff_vdw[:, :, :, 0], ff_vdw[:, :, :, 1], opt_dict["rotate"]) else: print("Loading Lennard-Jones force field from FFLJ_{x,y,z}") ff_vdw, lvec, _, atomic_info_or_head = io.load_vec_field("FFLJ", data_format=args.output_format) - ff_vdw[0, :, :, :], ff_vdw[1, :, :, :] = rotate_ff(ff_vdw[0, :, :, :], ff_vdw[1, :, :, :], opt_dict["rotate"]) + ff_vdw[:, :, :, 0], ff_vdw[:, :, :, 1] = rotate_ff(ff_vdw[:, :, :, 0], ff_vdw[:, :, :, 1], opt_dict["rotate"]) if charged_system: print("Loading electrostatic force field from FFel_{x,y,z}") ff_electrostatics, lvec, _, atomic_info_or_head = io.load_vec_field("FFel", data_format=args.output_format) - ff_electrostatics[0, :, :, :], ff_electrostatics[1, :, :, :] = rotate_ff(ff_electrostatics[0, :, :, :], ff_electrostatics[1, :, :, :], opt_dict["rotate"]) + ff_electrostatics[:, :, :, 0], ff_electrostatics[:, :, :, 1] = rotate_ff(ff_electrostatics[:, :, :, 0], ff_electrostatics[:, :, :, 1], opt_dict["rotate"]) if args.boltzmann or args.bI: print("Loading Boltzmann force field from FFboltz_{x,y,z}") ff_boltzman, lvec, _, atomic_info_or_head = io.load_vec_field("FFboltz", data_format=args.output_format) - ff_boltzman[0, :, :, :], ff_boltzman[1, :, :, :] = rotate_ff(ff_boltzman[0, :, :, :], ff_boltzman[1, :, :, :], opt_dict["rotate"]) + ff_boltzman[:, :, :, 0], ff_boltzman[:, :, :, 1] = rotate_ff(ff_boltzman[:, :, :, 0], ff_boltzman[:, :, :, 1], opt_dict["rotate"]) if applied_bias: print("Loading electrostatic contribution from applied bias from FFkpfm_t0sV_{x,y,z} and FFkpfm_tVs0_{x,y,z}") ff_kpfm_t0sv, lvec, _, atomic_info_or_head = io.load_vec_field("FFkpfm_t0sV", data_format=args.output_format) ff_kpfm_tvs0, lvec, _, atomic_info_or_head = io.load_vec_field("FFkpfm_tVs0", data_format=args.output_format) - ff_kpfm_t0sv[0, :, :, :], ff_kpfm_t0sv[1, :, :, :] = rotate_ff(ff_kpfm_t0sv[0, :, :, :], ff_kpfm_t0sv[1, :, :, :], opt_dict["rotate"]) - ff_kpfm_tvs0[0, :, :, :], ff_kpfm_tvs0[1, :, :, :] = rotate_ff(ff_kpfm_tvs0[0, :, :, :], ff_kpfm_tvs0[1, :, :, :], opt_dict["rotate"]) + ff_kpfm_t0sv[:, :, :, 0], ff_kpfm_t0sv[:, :, :, 1] = rotate_ff(ff_kpfm_t0sv[:, :, :, 0], ff_kpfm_t0sv[:, :, :, 1], opt_dict["rotate"]) + ff_kpfm_tvs0[:, :, :, 0], ff_kpfm_tvs0[:, :, :, 1] = rotate_ff(ff_kpfm_tvs0[:, :, :, 0], ff_kpfm_tvs0[:, :, :, 1], opt_dict["rotate"]) ff_kpfm_t0sv = ff_kpfm_t0sv * opt_dict["pol_s"] ff_kpfm_tvs0 = ff_kpfm_tvs0 * opt_dict["pol_t"] From bf93f8c342434957a3820f6c2065950a1a853409 Mon Sep 17 00:00:00 2001 From: Martin Ondracek Date: Wed, 28 May 2025 15:55:51 +0200 Subject: [PATCH 6/9] Rename RG to RGB in ppafm/PPPlot.py --- ppafm/PPPlot.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/ppafm/PPPlot.py b/ppafm/PPPlot.py index c042a7e1..07b01d91 100644 --- a/ppafm/PPPlot.py +++ b/ppafm/PPPlot.py @@ -47,7 +47,7 @@ def plotGeom(atoms=None, bonds=None, atomSize=default_atom_size): plotAtoms(atoms, atomSize=atomSize) -def colorize_XY2RG(Xs, Ys): +def colorize_XY2RGB(Xs, Ys): r = np.sqrt(Xs**2 + Ys**2) vmax = r[5:-5, 5:-5].max() Red = 0.5 * Xs / vmax + 0.5 @@ -120,7 +120,7 @@ def plotImages( plt.close() -def plotVecFieldRG( +def plotVecFieldRGB( prefix, dXs, dYs, @@ -148,7 +148,7 @@ def plotVecFieldRG( # print(" plotting ", i) write_plotting_slice(i) plt.figure(figsize=(10, 10)) - HSBs, vmax = colorize_XY2RG(dXs[i], dYs[i]) + HSBs, vmax = colorize_XY2RGB(dXs[i], dYs[i]) plt.imshow(HSBs, extent=extent, origin="lower", interpolation=interpolation) plotGeom(atoms, bonds, atomSize=atomSize) plt.xlabel(r" Tip_x $\AA$") From 8257dd9c293a8bb5404fc9f16d7f0ea612a4605c Mon Sep 17 00:00:00 2001 From: Martin Ondracek Date: Wed, 28 May 2025 16:19:05 +0200 Subject: [PATCH 7/9] Replace the erroneous expression ff_pauli[1, :, :, :, 1] by ff_pauli[:, :, :, 1] in ppafm/cli/relaxed_scan.py --- ppafm/cli/relaxed_scan.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ppafm/cli/relaxed_scan.py b/ppafm/cli/relaxed_scan.py index bd0b13a7..24907145 100755 --- a/ppafm/cli/relaxed_scan.py +++ b/ppafm/cli/relaxed_scan.py @@ -104,7 +104,7 @@ def main(argv=None): print("Loading Pauli force field from FFpauli_{x,y,z}") ff_pauli, lvec, _, atomic_info_or_head = io.load_vec_field("FFpauli", data_format=args.output_format) - ff_pauli[:, :, :, 0], ff_pauli[1, :, :, :, 1] = rotate_ff(ff_pauli[:, :, :, 0], ff_pauli[:, :, :, 1], opt_dict["rotate"]) + ff_pauli[:, :, :, 0], ff_pauli[:, :, :, 1] = rotate_ff(ff_pauli[:, :, :, 0], ff_pauli[:, :, :, 1], opt_dict["rotate"]) print("Loading vdW force field from FFvdW_{x,y,z}") ff_vdw, lvec, _, atomic_info_or_head = io.load_vec_field("FFvdW", data_format=args.output_format) From 24d488cb3c65fe4f13223adc86330bd1e0d4e30a Mon Sep 17 00:00:00 2001 From: Martin Ondracek Date: Thu, 29 May 2025 16:03:02 +0200 Subject: [PATCH 8/9] Use default "slices" argument to plot all slices with plotImages in ppafm/cli/plot_results.py --- ppafm/cli/plot_results.py | 10 ---------- 1 file changed, 10 deletions(-) diff --git a/ppafm/cli/plot_results.py b/ppafm/cli/plot_results.py index 4e79739c..b54162d3 100644 --- a/ppafm/cli/plot_results.py +++ b/ppafm/cli/plot_results.py @@ -141,7 +141,6 @@ def main(argv=None): dirname + "/xy" + atoms_str + cbar_str, pp_positions[:, :, :, 0], pp_positions[:, :, :, 1], - slices=list(range(0, pp_positions.shape[2])), BG=pp_positions[:, :, :, 2], extent=extent, atoms=atoms, @@ -164,7 +163,6 @@ def main(argv=None): PPPlot.plotImages( dirname + "/IETS" + atoms_str + cbar_str, iets, - slices=list(range(0, iets.shape[2])), zs=tip_positions_z, extent=extent, atoms=atoms, @@ -177,7 +175,6 @@ def main(argv=None): PPPlot.plotImages( dirname + "/Evib" + atoms_str + cbar_str, e_vib[:, :, :, 0], - slices=list(range(0, iets.shape[2])), zs=tip_positions_z, extent=extent, atoms=atoms, @@ -190,7 +187,6 @@ def main(argv=None): PPPlot.plotImages( dirname + "/Kvib" + atoms_str + cbar_str, 16.0217662 * eigenvalue_k[:, :, :, 0], - slices=(range(0, iets.shape[2])), zs=tip_positions_z, extent=extent, atoms=atoms, @@ -206,7 +202,6 @@ def main(argv=None): PPPlot.plotImages( dirname + "/OutI" + atoms_str + cbar_str, current, - slices=(range(0, current.shape[2])), zs=tip_positions_z, extent=extent, atoms=atoms, @@ -344,7 +339,6 @@ def main(argv=None): PPPlot.plotImages( dir_name_amplitude + "/df" + atoms_str + cbar_str, dfs, - slices=list(range(0, dfs.shape[2])), zs=tip_positions_z + parameters.Amplitude / 2.0, extent=extent, cmap=parameters.colorscale, @@ -371,7 +365,6 @@ def main(argv=None): PPPlot.plotImages( dir_name_amplitude + "/df_laplace" + atoms_str + cbar_str, df_laplace_filtered, - slices=list(range(0, len(dfs))), zs=tip_positions_z + parameters.Amplitude / 2.0, extent=extent, cmap=parameters.colorscale, @@ -405,7 +398,6 @@ def main(argv=None): PPPlot.plotImages( dir_name_lcpd + "/LCPD" + atoms_str + cbar_str, lcpd, - slices=list(range(0, lcpd.shape[2])), zs=tip_positions_z + parameters.Amplitude / 2.0, extent=extent, cmap=parameters.colorscale_kpfm, @@ -422,7 +414,6 @@ def main(argv=None): PPPlot.plotImages( dir_name_lcpd + "/_Asym-LCPD" + atoms_str + cbar_str, lcpd, - slices=list(range(0, lcpd.shape[2])), zs=tip_positions_z + parameters.Amplitude / 2.0, extent=extent, cmap=parameters.colorscale_kpfm, @@ -459,7 +450,6 @@ def main(argv=None): PPPlot.plotImages( dirname + "/Fz" + atoms_str + cbar_str, fzs, - slices=listlist(range(0, fzs.shape[2])), zs=tip_positions_z, # + parameters.Amplitude / 2.0, # no oscillation to the force extent=extent, cmap=parameters.colorscale, From 248047effa08c3ef5f238fac7e8df748f019e416 Mon Sep 17 00:00:00 2001 From: Martin Ondracek Date: Mon, 2 Jun 2025 15:21:32 +0200 Subject: [PATCH 9/9] Merge functions prepareArrays, setFF_shape, and setFF into a new setFF --- ppafm/HighLevel.py | 85 ++++++++++++++++++++++++------- ppafm/cli/generateTraining_PVE.py | 6 +-- ppafm/core.py | 50 +----------------- 3 files changed, 70 insertions(+), 71 deletions(-) diff --git a/ppafm/HighLevel.py b/ppafm/HighLevel.py index 92fd294a..fe310e9e 100755 --- a/ppafm/HighLevel.py +++ b/ppafm/HighLevel.py @@ -160,8 +160,7 @@ def perform_relaxation( FF += FFboltz if bFFtotDebug: io.save_vec_field("FFtotDebug", FF, lvec) - core.setFF_shape(np.shape(FF), lvec, parameters=parameters) - core.setFF_Fpointer(FF) + setFF(FF, lvec=lvec, parameters=parameters) if (np.array(parameters.stiffness) < 0.0).any(): parameters.stiffness = np.array([parameters.klat, parameters.klat, parameters.krad]) if verbose > 0: @@ -197,21 +196,71 @@ def perform_relaxation( # ==== Forcefield grid generation -def prepareArrays(FF, Vpot, parameters): - if parameters.gridN[2] <= 0: - PPU.autoGridN() - if FF is None: - gridN = parameters.gridN +def setFF(FF=None, computeVpot=False, n=None, lvec=None, parameters=None, verbose=True): + + # Find gridN + if FF is not None: + gridN = np.shape(FF)[0:3] + if parameters is not None: + parameters.gridN = gridN + else: + if n is not None: + gridN = np.array(n, dtype=np.int32) + elif parameters is not None: + if parameters.gridN[2] <= 0: + gridN = PPU.autoGridN(parameters) + else: + gridN = parameters.gridN + else: + raise ValueError("FF dimensions not set !!") + # Create a new array for FF if needed FF = np.zeros((gridN[0], gridN[1], gridN[2], 3)) + if verbose: + print("setFF() gridN: ", gridN) + core.setGridN(gridN) + + # Set pointer to FF + if len(FF.shape) == 4 and FF.shape[-1] == 3: + if verbose: + print("setFF() Creating a pointer to a vector field") + core.setFF_Fpointer(FF) + elif len(FF.shape) == 3 or FF.shape[-1] == 1: + if verbose: + print("setFF() Creating a pointer to a scalar field") + core.setFF_Epointer(FF) else: - gridN = np.shape(FF) - parameters.gridN = gridN - core.setFF_Fpointer(FF) - if Vpot: + raise ValueError("setFF: Array dimensions wrong for both vector and array field !!") + + # Create a scalar (potential) field if required + if computeVpot: + if verbose: + print("setFF() Creating a pointer to a scalar field") V = np.zeros((gridN[0], gridN[1], gridN[2])) core.setFF_Epointer(V) else: V = None + + # Set lattice vectors or grid geometry + if (lvec is None) and (parameters is not None): + lvec = np.array( + [ + parameters.gridA, + parameters.gridB, + parameters.gridC, + ], + dtype=np.float64, + ).copy() + lvec = np.array(lvec, dtype=np.float64) + if lvec.shape == (3, 3): + lvec = lvec.copy() + elif lvec.shape == (4, 3): + lvec = lvec[1:, :].copy() + else: + raise ValueError("lvec matrix has a wrong format !!") + if verbose: + print("setFF() lvec: ", lvec) + core.setGridCell(lvec) + return FF, V @@ -231,10 +280,9 @@ def computeLJ(geomFile, speciesFile, geometry_format=None, save_format=None, com # --- prepare LJ parameters iPP = PPU.atom2iZ(parameters.probeType, elem_dict) # --- prepare arrays and compute - FF, V = prepareArrays(None, computeVpot, parameters=parameters) + FF, V = setFF(None, computeVpot, lvec=lvec, parameters=parameters) if verbose > 0: print("FFLJ.shape", FF.shape) - core.setFF_shape(np.shape(FF), lvec, parameters=parameters) # shift atoms to the coordinate system in which the grid origin is zero Rs0 = shift_positions(Rs, -lvec[0]) @@ -304,8 +352,7 @@ def computeDFTD3(input_file, df_params="PBE", geometry_format=None, save_format= coeffs = core.computeD3Coeffs(Rs, iZs, iPP, df_params) # Compute the force field - FF, V = prepareArrays(None, compute_energy, parameters=parameters) - core.setFF_shape(np.shape(FF), lvec, parameters=parameters) + FF, V = setFF(None, compute_energy, lvec=lvec, parameters=parameters) core.getDFTD3FF(shift_positions(Rs, -lvec[0]), coeffs) # Save to file @@ -336,8 +383,7 @@ def computeELFF_pointCharge(geomFile, geometry_format=None, tip="s", save_format if verbose > 0: print(parameters.gridN, parameters.gridA, parameters.gridB, parameters.gridC) _, Rs, Qs = PPU.parseAtoms(atoms, elem_dict=elem_dict, autogeom=False, PBC=parameters.PBC, lvec=lvec, parameters=parameters) - FF, V = prepareArrays(None, computeVpot, parameters=parameters) - core.setFF_shape(np.shape(FF), lvec, parameters=parameters) + FF, V = setFF(None, computeVpot, lvec=lvec, parameters=parameters) # shift atoms to the coordinate system in which the grid origin is zero Rs0 = shift_positions(Rs, -lvec[0]) @@ -446,8 +492,9 @@ def subtractCoreDensities( print("V : ", V, " N: ", N, " dV: ", dV) if verbose > 0: print("sum(RHO): ", rho.sum(), " Nelec: ", rho.sum() * dV, " voxel volume: ", dV) # check sum - core.setFF_shape(rho.shape, lvec, parameters=parameters) # set grid sampling dimension and shape - core.setFF_Epointer(rho) # set pointer to array with density data (to write into) + + # set grid sampling dimension and shape + setFF(rho, lvec=lvec, parameters=parameters) if verbose > 0: print(">>> Projecting Core Densities ... ") diff --git a/ppafm/cli/generateTraining_PVE.py b/ppafm/cli/generateTraining_PVE.py index 7e897a04..7782878e 100644 --- a/ppafm/cli/generateTraining_PVE.py +++ b/ppafm/cli/generateTraining_PVE.py @@ -9,7 +9,7 @@ from ppafm import io from .. import common, core -from ..HighLevel import prepareArrays, relaxedScan3D +from ..HighLevel import relaxedScan3D, setFF file_format = "xsf" @@ -34,10 +34,8 @@ parameters.gridC = lvec_t[3] # must be before parseAtoms print(parameters.gridN, parameters.gridA, parameters.gridB, parameters.gridC) -force_field, _ = prepareArrays(None, False) - +force_field, _ = setFF(None, False, lvec=lvec_t, parameters=parameters) print("FFLJ.shape", force_field.shape) -core.setFF_shape(np.shape(force_field), lvec_t, parameters=parameters) base_dir = os.getcwd() paths = ["out1", "out2"] diff --git a/ppafm/core.py b/ppafm/core.py index b84276b0..929778e1 100755 --- a/ppafm/core.py +++ b/ppafm/core.py @@ -42,34 +42,11 @@ def setGridN(n): lib.setGridCell.restype = None -def setGridCell(cell=None, parameters=None): - if cell is None: - cell = np.array( - [ - parameters.gridA, - parameters.gridB, - parameters.gridC, - ], - dtype=np.float64, - ).copy() - cell = np.array(cell, dtype=np.float64) - if cell.shape == (3, 3): - cell = cell.copy() - elif cell.shape == (4, 3): - cell = cell[1:, :].copy() - else: - raise ValueError("cell has wrong format") - exit() - print("cell", cell) +def setGridCell(cell): + cell = np.array(cell[-3:, 0:3], copy=True, dtype=np.float64) lib.setGridCell(cell) -def setFF_shape(n_, cell, parameters=None): - n = np.array(n_).astype(np.int32) - lib.setGridN(n) - setGridCell(cell, parameters=parameters) - - # void setFF_pointer( double * gridF, double * gridE ) lib.setFF_Fpointer.argtypes = [array4d] lib.setFF_Fpointer.restype = None @@ -108,29 +85,6 @@ def deleteFF_Epointer(): lib.deleteFF_Epointer() -def setFF(cell=None, gridF=None, gridE=None, parameters=None): - n_ = None - if gridF is not None: - setFF_Fpointer(gridF) - n_ = np.shape(gridF) - if gridE is not None: - setFF_Epointer(gridE) - n_ = np.shape(gridF) - if cell is None: - cell = np.array( - [ - parameters.gridA, - parameters.gridB, - parameters.gridC, - ] - ).copy() - if n_ is not None: - print("setFF() n_ : ", n_) - setFF_shape(n_, cell, parameters=parameters) - else: - "Warrning : setFF shape not set !!!" - - # void setRelax( int maxIters, double convF2, double dt, double damping ) lib.setRelax.argtypes = [c_int, c_double, c_double, c_double] lib.setRelax.restype = None