ppafm.ocl

ppafm.ocl.AFMulator

class ppafm.ocl.AFMulator.AFMulator(pixPerAngstrome=10, lvec=None, scan_dim=(128, 128, 30), scan_window=((2.0, 2.0, 7.0), (18.0, 18.0, 10.0)), iZPP=8, QZs=[0.1, 0, -0.1, 0], Qs=[-10, 20, -10, 0], rho=None, sigma=0.71, rho_delta=None, A_pauli=18.0, B_pauli=1.0, fdbm_vdw_type='D3', d3_params='PBE', lj_vdw_damp=2, df_steps=10, tipR0=[0.0, 0.0, 3.0], tipStiffness=[0.25, 0.25, 0.0, 30.0], npbc=(1, 1, 0), f0Cantilever=30300, kCantilever=1800, colorscale='gray', minimize_memory=False)

Bases: object

Simulate Atomic force microscope images of molecules.

Parameters:

• pixPerAngstrome – int. Number of pixels (voxels) per angstrom in force field grid.

• lvec – np.ndarray of shape (4, 3) or None. Unit cell boundaries for force field. First (row) vector specifies the origin, and the remaining three vectors specify the edge vectors of the unit cell. If None, will be calculated automatically from scan_window and tipR0, leaving some additional space on each side.

• scan_dim – tuple of three ints. Number of scan steps in the (x, y, z) dimensions. The size of the resulting images have a size of scan_dim[0] x scan_dim[1] and the number images at different heights is scan_dim[2] - df_steps + 1.

• scan_window – tuple ((x_min, y_min, z_min), (x_max, y_max, z_max)). The minimum and maximum coordinates of scan region in angstroms. The scan region is a rectangular box with the opposing corners defined by the coordinates in scan_window. Note that the step size in z dimension is the scan window size in z divided by scan_dim[2], and the scan in z-direction proceeds for scan_dim[2] steps, so the final step is one step short of z_min.

• iZPPs – int. Element of probe particle.

• QZs – Array of length 4. Positions of tip charges along z-axis relative to probe-particle center in angstroms.

• Qs – Array of length 4. Values of tip charges in units of e. Some can be zero.

• rho – Dict or TipDensity. Tip charge density. Used with Hartree potentials. Overrides QZs and Qs. The dict should contain float entries for at least of one the following ‘s’, ‘px’, ‘py’, ‘pz’, ‘dz2’, ‘dy2’, ‘dx2’, ‘dxy’, ‘dxz’, ‘dyz’. The tip charge density will be a linear combination of the specified multipole types with the specified weights.

• sigma – float. Width of tip density distribution if rho is a dict of multipole coefficients.

• rho_delta – TipDensity. Tip delta charge density. Used in FDBM approximation for calculating the electrostatic interaction force.

• A_pauli – float. Prefactor for Pauli repulsion when using the FDBM.

• B_pauli – float. Exponent for Pauli repulsion when using the FDBM.

• fdbm_vdw_type – ‘D3’ or ‘LJ’. Type of vdW interaction to use with the FDBM. ‘D3’ is for Grimme-D3 and ‘LJ’ uses standard Lennard-Jones vdW.

• d3_params – str or dict. Functional-specific scaling parameters for Grimme-D3. Can be a str with the functional name or a dict with manually specified parameters. Used in FDBM. See add_dftd3() for further explanation.

• lj_vdw_damp – int. Type of damping to use in vdw calculation for FDBM when fdbm_vdw_type==’LJ’. -1: no damping, 0: constant, 1: R2, 2: R4, 3: invR4, 4: invR8.

• df_steps – int. Number of steps in z convolution. The total amplitude is df_steps times scan z-step size.

• tipR0 – array of length 3. Probe particle equilibrium position (x, y, z) in angstroms.

• tipStiffness – array of length 4. Harmonic spring constants (x, y, z, r) in N/m for holding the probe particle to the tip.

• npbc – tuple of three ints. How many periodic images of atoms to use in (x, y, z) dimensions. Used for calculating Lennard-Jones force field and electrostatic field from point charges. Electrostatic field from a Hartree potential defined on a grid is always considered to be periodic.

• f0Cantilever – float. Resonance frequency of the cantilever in Hz.

• kCantilever – float. Harmonic spring constant of the cantilever in N/m.

• colorscale – str. Colorscale for output images.

• minimize_memory – bool. Release device memory as soon as it’s not needed. Can help save some memory, but can also make the simulation significantly slower when run in a loop where parameters change between iterations.

bMergeConv = False

bSaveFF = False

check_scan_window()

Check that scan window does not extend beyond any non-periodic boundaries.

eval(xyzs, Zs, qs, rho_sample=None, sample_lvec=None, rot=array([[1., 0., 0.], [0., 1., 0.], [0., 0., 1.]]), rot_center=None, REAs=None, X=None, plot_to_dir=None)

Prepare and evaluate AFM image.

Parameters:

• xyzs – np.ndarray of shape (num_atoms, 3). Positions of atoms in x, y, and z.

• Zs – np.ndarray of shape (num_atoms,). Elements of atoms.

• qs – np.ndarray of shape (num_atoms,) or HartreePotential or None. Charges of atoms or hartree potential. If None, then no electrostatics are used.

• rho_sample – ElectronDensity or None. Sample electron density. If not None, then FDBM is used for calculating the Pauli repulsion force.

• sample_lvec – np.ndarray of shape (3, 3) or None. Unit cell lattice vectors for periodic images of atoms. If None, periodic boundaries are disabled, unless qs is HartreePotential and the lvec from the Hartree potential is used instead. If npbc = (0, 0, 0), then has no function.

• rot – np.ndarray of shape (3, 3). Rotation matrix to apply to atom positions.

• rot_center – np.ndarray of shape (3,). Center for rotation. Defaults to center of atom coordinates.

• REAs – np.ndarray of shape (num_atoms, 4). Lennard Jones interaction parameters. Calculated automatically if None.

• X – np.ndarray of shape (self.scan_dim[0], self.scan_dim[1], self.scan_dim[2]-self.df_steps+1)). Array where AFM image will be saved. If None, will be created automatically.

• plot_to_dir – str or None. If not None, plot the generated AFM images to this directory.

Returns:

X – np.ndarray. Output AFM images. If input X is not None, this is the same array object as X with values overwritten.

evalAFM(X=None)

Evaluate AFM image. Run after preparing force field and scanner.

Parameters:

X – np.ndarray of shape (self.scan_dim[0], self.scan_dim[1], self.scan_dim[2]-self.df_steps+1)). Array where AFM image will be saved. If None, will be created automatically.

Returns:

X – np.ndarray. Output AFM images. If input X is not None, this is the same array object as X with values overwritten.

eval_(mol)

classmethod from_params(file_path='./params.ini')

Construct an AFMulator instance from a params.ini file.

Parameters:

file_path – str. Path to the params.ini file to load.

load_params(file_path='./params.ini')

Update the parameters of this AFMulator instance from a params.ini file.

Parameters:

file_path – str. Path to the params.ini file to load.

plot_images(X, outdir='afm_images', prefix='df', **plot_kwargs)

Plot simulated AFM images and save them to a directory.

Parameters:

• X – np.ndarray. AFM images to plot.

• outdir – str. Path to directory where files are saved.

• prefix – str. Prefix string for saved files.

• plot_kwargs – Additional keyword arguments passed to plotImages().

prepareFF(xyzs, Zs, qs, rho_sample=None, sample_lvec=None, rot=array([[1., 0., 0.], [0., 1., 0.], [0., 0., 1.]]), rot_center=None, REAs=None)

Prepare molecule parameters and calculate force field.

Parameters:

• xyzs – np.ndarray of shape (num_atoms, 3). Positions of atoms in x, y, and z.

• Zs – np.ndarray of shape (num_atoms,). Elements of atoms.

• qs – np.ndarray of shape (num_atoms,) or HartreePotential or None. Charges of atoms or hartree potential. If None, then no electrostatics are used.

• rho_sample – ElectronDensity or None. Sample electron density. If not None, then FDBM is used for calculating the Pauli repulsion force. Requires rho_delta to be set and qs has to be HartreePotential.

• sample_lvec – np.ndarray of shape (3, 3) or None. Unit cell lattice vectors for periodic images of atoms. If None, periodic boundaries are disabled, unless qs is HartreePotential and the lvec from the Hartree potential is used instead. If npbc = (0, 0, 0), then has no function.

• rot – np.ndarray of shape (3, 3). Rotation matrix to apply to atom positions.

• rot_center – np.ndarray of shape (3,). Center for rotation. Defaults to center of atom coordinates.

• REAs – np.ndarray of shape (num_atoms, 4). Lennard Jones interaction parameters. Calculated automatically if None.

prepareScanner()

Prepare scanner. Run after preparing force field.

relaxParams = [0.5, 0.1, 0.020000000000000004, 0.5]

saveDebugXSF_FF(fname, F)

saveFF()

save_params(file_path='./params.ini')

Save the parameters of this AFMulator instance to a params.ini file.

Parameters:

file_path – str. Path to the file where parameters are saved.

setBPauli(B_pauli=1.0)

Set Pauli repulsion exponent used in FDBM.

setLvec(lvec=None, pixPerAngstrome=None)

Set forcefield lattice vectors. If lvec is not given it is inferred from the scan window.

setQs(Qs, QZs)

Set tip point charges.

setRho(rho=None, sigma=0.71, B_pauli=1.0)

Set tip charge distribution.

Parameters:

• rho – Dict, TipDensity, or None. Tip charge density. If None, the existing density is deleted.

• sigma – float. Tip charge density distribution when rho is a dict.

• B_pauli – float. Pauli repulsion exponent for tip density when using FDBM.

setRhoDelta(rho_delta=None)

Set tip electron delta-density that is used for electrostatic interaction in FDBM.

Parameters:

rho_delta – TipDensity or None. Tip electron delta-density. If None, the existing density is deleted.

setScanWindow(scan_window=None, scan_dim=None, df_steps=None)

Set scanner scan window.

setStiffness(tipStiffness)

Set the tip bending stiffness.

ppafm.ocl.AFMulator.get_lvec(scan_window, pad=(3.0, 3.0, 5.0), tipR0=(0.0, 0.0, 3.0), pixPerAngstrome=10)

ppafm.ocl.AFMulator.quick_afm(file_path, scan_size=(16, 16), offset=(0, 0), distance=8.0, scan_step=(0.1, 0.1, 0.1), probe_type=8, charge=-0.1, tip='dz2', sigma=0.71, num_heights=10, amplitude=1.0, out_dir=None)

Make an AFM simulation from a .cube, .xsf, or .xyz file, and print images to a folder.

Parameters:

• file_path – str. Path to input file.

• scan_size – tuple of length 2. Size of scan region in angstroms.

• offset – tuple of length 2. Offset to center of scan region in angstroms.

• distance – float. Furthest distance of probe tip from sample.

• scan_step – tuple of length 3. Scan steps in x, y, and z dimensions.

• probe_type – int. Probe particle type.

• charge – float. Tip charge.

• tip – str. Tip multipole type.

• sigma – float. Width of tip charge distribution.

• num_heights – int. Number of different heights to scan.

• amplitude – float. Oscillation amplitude in angstroms.

• out_dir – str or None. Output folder path. If None, defaults to “./afm_” + input_file_name.

ppafm.ocl.field

class ppafm.ocl.field.AtomProjection

Bases: object

to generate reference output maps ( Ys ) in generator for Neural Network training

Rfunc = None

Rmax = 10.0

Rpp = 2.0

drStep = 0.1

dzmax = 2.0

dzmax_s = inf

elipticity = 0.5

makeCoefsZR(Zs, ELEMENTS)

make atomic coeficients used e.g. in MultiMap

prepareBuffers(atoms, prj_dim, coefs=None, bonds2atoms=None, Rfunc=None, elem_channels=None)

allocate GPU buffers

releaseBuffers()

deallocated all GPU buffers

run_evalAtomRfunc(poss=None, Eout=None, tipRot=None, local_size=(32,))

kernel to produce multiple channels of vdW Sphere maps each coresponding to different atom type

run_evalBondEllipses(poss=None, Eout=None, tipRot=None, local_size=(32,))

kernel to produce multiple channels of vdW Sphere maps each coresponding to different atom type

run_evalCoulomb(poss=None, Eout=None, local_size=(32,))

kernel producing coulomb potential and field

run_evalHartreeGradient(pot, poss=None, Eout=None, h=None, rot=array([[1., 0., 0.], [0., 1., 0.], [0., 0., 1.]]), rot_center=None, local_size=(32,))

Get electric field as the negative gradient of a Hartree potential.

Parameters:

• pot – HartreePotential. Hartree potential to differentiate.

• poss – np.ndarray or None. Position grid for points to get the field at.

• Eout – np.ndarray or None. Output array. If None, will be created automatically.

• h – float > 0.0 or None. Finite difference step size (one-sided) in angstroms. If None, the default value DEFAULT_FD_STEP is used.

• rot – np.ndarray of shape (3, 3). Rotation matrix to apply to the position coordinates.

• rot_center – np.ndarray of shape (3,). Point around which rotation is performed.

• local_size – tuple of a single int. Size of local work group on device.

run_evalLorenz(poss=None, Eout=None, local_size=(32,))

kernel producing lorenzian function around each atom

run_evalMultiMapSpheres(poss=None, Eout=None, tipRot=None, bOccl=0, Rmin=1.4, Rstep=0.1, local_size=(32,))

kernel to produce multiple channels of vdW Sphere maps each containing atoms with different vdW radius

run_evalMultiMapSpheresElements(poss=None, Eout=None, tipRot=None, bOccl=0, local_size=(32,))

kernel to produce multiple channels of vdW Sphere maps each containing atoms with different vdW radius

run_evalQdisks(poss=None, Eout=None, tipRot=None, local_size=(32,))

kernel producing atoms disks with positive and negative value encoding charge

run_evalSphereCaps(poss=None, Eout=None, tipRot=None, local_size=(32,))

kernel producing spherical caps (just to top most part of vdW sphere)

run_evalSpheres(poss=None, Eout=None, tipRot=None, local_size=(32,))

kernel producing van der Waals spheres

run_evalSpheresType(poss=None, Eout=None, tipRot=None, bOccl=0, local_size=(32,))

kernel to produce multiple channels of vdW Sphere maps each coresponding to different atom type

run_evaldisks(poss=None, Eout=None, tipRot=None, offset=0.0, local_size=(32,))

kernel producing atomic disks with conical profile

run_evaldisks_occlusion(poss=None, Eout=None, tipRot=None, local_size=(32,))

kernel producing atomic disks occluded by higher nearby atoms

setAtomTypes(types, sel=[1, 6, 8])

setup selection of atomic types for SpheresType kernel and upload them to GPU

tgMax = 0.5

tgWidth = 0.1

tryReleaseBuffers()

deallocated all GPU buffers (those which exists)

updateBuffers(atoms=None, coefs=None, poss=None)

upload data to GPU

zmargin = 0.2

zmin = -3.0

class ppafm.ocl.field.D3Params(ctx)

Bases: object

pyopencl device buffer handles to Grimme-D3 parameters. Each buffer is allocated on first access.

Parameters:

ctx – pyopencl.Context. OpenCL context for device buffer. Defaults to oclu.ctx.

property cl_r4r2

See R4R2.

property cl_rcov

See R_COV.

property cl_rcut

See load_R0().

property cl_ref_c6

See load_ref_c6().

property cl_ref_cn

See REF_CN.

release()

Release device buffers.

class ppafm.ocl.field.DataGrid(array, lvec, shape=None, ctx=None)

Bases: object

Class for holding data on a grid. The data can be stored on the CPU host or an OpenCL device.

Parameters:

• array – np.ndarray or pyopencl.Buffer. Array values on a 3D grid with possibly multiple components.

• lvec – array-like of shape (4, 3). Unit cell boundaries. First (row) vector specifies the origin, and the remaining three vectors specify the edge vectors of the unit cell.

• shape – array-like of length 3 or 4. Grid shape when array is a pyopencl.Buffer.

• ctx – pyopencl.Context. OpenCL context for device buffer. Defaults to oclu.ctx.

add_mult(array, scale=1.0, in_place=True, local_size=(32,), queue=None)

Multiply the values of another data grid and add them to the values of this data grid.

Parameters:

• array – DataGrid. Grid whose values to scale and add.

• scale – float. Value by which values in array are multiplied.

• in_place – bool. Whether to do operation in place or to create a new array.

• local_size – tuple of a single int. Size of local work group on device.

• queue – pyopencl.CommandQueue. OpenCL queue on which operation is performed. Defaults to oclu.queue.

Returns:

grid_out – Same type as self. New data grid with result.

property array

Host array as np.ndarray. If the grid currently only exists on the device, it is copied to the host memory.

property cell_vol

The volume of a grid cell in angstrom^3.

property cl_array

Device array as pyopencl.buffer. If the grid currently only exists on the host, it is copied to the device memory.

clamp(minimum=-inf, maximum=inf, clamp_type='hard', soft_clamp_width=1.0, in_place=True, local_size=(32,), queue=None)

Clamp data grid values to a specified range. The 'hard' clamp simply clips values that are out of range, and the 'soft' clamp uses a sigmoid to smoothen the transition.

Parameters:

• minimum – float. Values below minimum are set to minimum.

• maximum – float. Values above maximum are set to maximum.

• clamp_type – str. Type of clamp to use: 'soft' or 'hard'.

• soft_clamp_width – float. Width of transition region for soft clamp.

• in_place – bool. Whether to do operation in place or to create a new array.

• local_size – tuple of a single int. Size of local work group on device.

• queue – pyopencl.CommandQueue. OpenCL queue on which operation is performed. Defaults to oclu.queue.

Returns:

grid_out – Same type as self. New data grid with result.

classmethod from_file(file_path, scale=1.0)

Load grid data and atoms from a .cube or .xsf file.

Parameters:

• file_path – str. Path to file to load.

• scale – float. Scaling factor for the returned data grid values.

Returns:

• data – class type. Data grid object.

• xyzs – np.ndarray of shape (num_atoms, 3). Atom coordinates.

• Zs – np.ndarray of shape (num_atoms,). Atomic numbers.

grad(scale=1.0, array_out=None, order='C', local_size=(32,), queue=None)

Get the centered finite difference gradient of the data grid. Uses periodic boundary conditions at the edge of the grid.

The datagrid has to be either 3D, or 4D with self.shape[3] == 1.

The resulting array adds a 4th dimension with size 4 to the grid such that at indices 0-2 are the partial derivatives in x, y, and z directions, respectively, and at index 3 is the original scalar field.

Parameters:

• scale – float or array-like of size 4. Additional scaling factor for the output.

• array_out – pyopencl.Buffer or None. Output array. If None, then is created automatically.

• order – str, ‘C’ or ‘F’. Whether to save values in C or Fortran order.

• local_size – tuple of a single int. Size of local work group on device.

• queue – pyopencl.CommandQueue. OpenCL queue on which operation is performed. Defaults to oclu.queue.

Returns:

grid_out – DataGrid. New data grid with result.

interp_at(lvec_new, shape_new, array_out=None, rot=array([[1., 0., 0.], [0., 1., 0.], [0., 0., 1.]]), rot_center=array([0., 0., 0.]), local_size=(32,), queue=None)

Interpolate grid values onto another grid. Uses periodic boundary conditions.

Parameters:

• lvec_new – array-like of shape (4, 3). Unit cell boundaries for new grid.

• shape_new – array-like of length 3. New grid shape.

• array_out – pyopencl.Buffer or None. Output array. If None, then is created automatically.

• rot – np.ndarray of shape (3, 3). Rotation matrix to apply.

• rot_center – np.ndarray of shape (3,). Point around which rotation is performed.

• local_size – tuple of a single int. Size of local work group on device.

• queue – pyopencl.CommandQueue. OpenCL queue on which operation is performed. Defaults to oclu.queue.

Returns:

grid_out – Same type as self. New data grid with result.

power_positive(p=1.2, normalize=True, in_place=True, local_size=(32,), queue=None)

Raise every positive element in the grid into a power. Negative values are set to zero.

Parameters:

• p – float. Power to rise to.

• normalize – bool. Whether to normalize the values after setting negative values to zero. The normalization is done by scaling the values such that the total sum of the values in the array remains unchanged after eliminating the negative values.

• in_place – bool. Whether to do operation in place or to create a new array.

• local_size – tuple of a single int. Size of local work group on device.

• queue – pyopencl.CommandQueue. OpenCL queue on which operation is performed. Defaults to oclu.queue.

Returns:

grid_out – Same type as self. New data grid with result.

release(keep_on_host=True)

Release device buffer.

Arguments: keep_on_host: bool. If the grid currently only exists on the device, it is copied to the host memory before release.

property step

Array of vectors pointing single steps along the grid for each lattice vector.

to_file(file_path, clamp=None)

Save data grid to file(s).

Supported file types are .xsf and .npy.

Parameters:

• file_path – str. Path to saved file. For a 4D data grid, letters x, y, z, w are appended to the file path for each component, respectively.

• clamp – float or None. If not None, all values greater than this are clamped to this value.

update_array(array, lvec)

Update array contents. If the new array is the same size or smaller than the current array, the data is updated without a reallocation on the device.

class ppafm.ocl.field.ElectronDensity(array, lvec, shape=None, ctx=None)

Bases: DataGrid

Sample electron density. Units should be in e/Å^3.

class ppafm.ocl.field.FFTCrossCorrelation(rho, queue=None)

Bases: object

Do circular cross-correlation of sample Hartree potential or electron density with tip charge density via FFT.

Parameters:

• rho – TipDensity. Tip charge density.

• queue – pyopencl.CommandQueue. OpenCL queue on which operations are performed. Defaults to oclu.queue.

correlate(array, E=None, scale=1)

Cross-correlate input array with tip charge density.

Parameters:

• array – DataGrid or pyopencl.Buffer. Sample potential/density to cross-correlate with tip density. Has to be the same shape as rho.

• E – DataGrid or None. Output data grid. If None, is created automatically. The automatically created datagrid has the same lvec as self.rho.lvec.

• scale – float. Additional scaling factor for the output.

Returns:

E – DataGrid. Result of cross-correlation.

class ppafm.ocl.field.ForceField_LJC

Bases: object

Evaluate Lennard-Jones based force fields on an OpenCL device.

addLJ(local_size=(32,))

Add Lennard-Jones force and energy to the current force field grid.

Parameters:

local_size – tuple of a single int. Size of local work group on device.

add_dftd3(params='PBE', local_size=(64,))

Add van der Waals force and energy to the force field grid using the DFT-D3 method. Uses the Becke-Johnson damping method. Mainly useful in conjunction with the full-density based model.

The DFT-D3 parameters are adjusted based on the DFT functional. There are predefined scaling parameters for the following functionals: PBE, B1B95, B2GPPLYP, B3PW91, BHLYP, BMK, BOP, BPBE, CAMB3LYP, LCwPBE, MPW1B95, MPWB1K, mPWLYP, OLYP, OPBE, oTPSS, PBE38, PBEsol, PTPSS, PWB6K, revSSB, SSB, TPSSh, HCTH120, B2PLYP, B3LYP, B97D, BLYP, BP86, DSDBLYP, PBE0, PBE, PW6B95, PWPB95, revPBE0, revPBE38, revPBE, rPW86PBE, TPSS0, TPSS. See also DF_DEFAULT_PARAMS.

Otherwise, the parameters can be manually specified in a dict with the following entries:

• ‘s6’: Scaling parameter for r^-6 term.

• ‘s8’: Scaling parameter for r^-8 term.

• ‘a1’: Scaling parameter for cutoff radius.

• ‘a2’: Additive parameter for cutoff radius. Unit should be Bohr.

Parameters:

• params – str or dict. Functional-specific scaling parameters. Can be a str with the functional name or a dict with manually specified parameters.

• local_size – tuple of a single int. Size of local work group on device.

addvdW(damp_method=0, local_size=(32,))

Add Lennard-Jones van der Waals force and energy to the current force field grid.

Parameters:

• damp_method – int. Type of damping to use. -1: no damping, 0: constant, 1: R2, 2: R4, 3: invR4, 4: invR8.

• local_size – tuple of a single int. Size of local work group on device.

calc_force_fdbm(A=18.0, B=1.0, vdw_type='D3', d3_params='PBE', lj_vdw_damp=2, FE=None, rot=array([[1., 0., 0.], [0., 1., 0.], [0., 0., 1.]]), rot_center=array([0., 0., 0.]), local_size=(32,), bCopy=True, bFinish=True)

Calculate force field using the full density-based model.

Parameters:

• A – float. Prefactor for Pauli repulsion.

• B – float. Exponent used for Pauli repulsion.

• vdw_type – 'D3' or 'LJ'. Type of vdW interaction to use with the FDBM. 'D3' is for Grimme-D3 and 'LJ' uses standard Lennard-Jones vdW.

• d3_params – str or dict. Functional-specific scaling parameters for DFT-D3. Can be a str with the functional name or a dict with manually specified parameters. See add_dftd3().

• lj_vdw_damp – int. Type of damping to use in vdw calculation fdbm_vdw_type=='LJ'. -1: no damping, 0: constant, 1: R2, 2: R4, 3: invR4, 4: invR8.

• FE – np.ndarray or None. Array where output force field is copied to if bCopy==True. If None and bCopy==True, will be created automatically.

• rot – np.ndarray of shape (3, 3). Rotation matrix applied to the atom coordinates.

• rot_center – np.ndarray of shape (3,). Point around which rotation is performed.

• local_size – tuple of a single int. Size of local work group on device.

• bCopy – Bool. Whether to copy the calculated forcefield field to host.

• bFinish – Bool. Whether to wait for execution to finish.

Returns:

FE – np.ndarray if bCopy==True or None otherwise. Calculated force field and energy.

calc_force_hartree(FE=None, rot=array([[1., 0., 0.], [0., 1., 0.], [0., 0., 1.]]), rot_center=array([0., 0., 0.]), local_size=(32,), bCopy=True, bFinish=True)

Calculate force field for LJ + Hartree cross-correlated with tip density via FFT.

Parameters:

• FE – np.ndarray or None. Array where output force field is copied to if bCopy == True. If None and bCopy == True, will be created automatically.

• rot – np.ndarray of shape (3, 3). Rotation matrix applied to the atom coordinates.

• rot_center – np.ndarray of shape (3,). Point around which rotation is performed.

• local_size – tuple of a single int. Size of local work group on device.

• bCopy – Bool. Whether to copy the calculated forcefield field to host.

• bFinish – Bool. Whether to wait for execution to finish.

Returns:

FE – np.ndarray if bCopy==True or None otherwise. Calculated force field and energy.

downloadFF(FE=None)

Get the force field array from the device.

Parameters:

FE – np.ndarray or None. Array where output force field is copied to. If None, will be created automatically.

Returns:

FE – np.ndarray. Force field and energy.

initPoss(poss=None, nDim=None, lvec=None, pixPerAngstrome=10)

initSampling(lvec, pixPerAngstrome=10, nDim=None)

initialize(value=0, bFinish=False)

Initialize the force field to a constant value.

Parameters:

• value – float. Value assigned to every element of the force field grid.

• bFinish – Bool. Whether to wait for execution to finish.

makeFF(xyzs, cLJs, REAs=None, Zs=None, method='point-charge', FE=None, qs=None, pot=None, rho_sample=None, rho=None, rho_delta=None, A=18.0, B=1.0, fdbm_vdw_type='D3', d3_params='PBE', lj_vdw_damp=2, rot=array([[1., 0., 0.], [0., 1., 0.], [0., 0., 1.]]), rot_center=array([0., 0., 0.]), local_size=(32,), bRelease=True, bCopy=True, bFinish=True)

Generate the force field for a tip-sample interaction.

There are several methods for generating the force field:

• 'point-charge': Lennard-Jones + point-charge electrostatics for both tip and sample.

• 'hartree': Lennard-Jones + sample hartree potential cross-correlated with tip charge density for electrostatic interaction.

• 'fdbm': Approximated full density-based model. Pauli repulsion is calculated by tip-sample electron density overlap + attractive vdW like in Lennard-Jones. Electrostatic interaction is same as in ‘hartree’, except tip delta-density is used instead.

If pot, rho, or rho_delta is None and is required for the specified method, it has to be initialized beforehand with prepareBuffers().

Parameters:

• xyzs – np.ndarray of shape (n_atoms, 3). xyz positions.

• cLJs – np.ndarray of shape (n_atoms, 2). Lennard-Jones interaction parameters in AB form for each atom.

• REAs – np.ndarray of shape (n_atoms, 4) or None. Lennard-Jones interaction parameters in RE form for each atom. Required when method is ‘fdbm’, fdbm_vdw_type is ‘LJ’, and vdw_damp_method >= 1.

• Zs – np.ndarray of shape (n_atoms,). Atomic numbers. Required when method is ‘fdbm’.

• method – ‘point-charge’, ‘hartree’ or ‘fdbm’. Method for generating the force field.

• FE – np.ndarray or None. Array where output force field is copied to if bCopy == True. If None and bCopy == True, will be created automatically.

• qs – np.ndarray of shape (n_atoms,) or None. Point charges of atoms. Used when method is ‘point-charge’.

• pot – HartreePotential or None. Hartree potential used for electrostatic interaction when method is ‘hartree’ or ‘fdbm’.

• rho_sample – ElectronDensity or None. Sample electron density. Used for Pauli repulsion when method is ‘fdbm’.

• rho – TipDensity or None. Probe tip charge density. Used for electrostatic interaction when method is ‘hartree’ and Pauli repulsion when method is ‘fdbm’.

• rho_delta – TipDensity or None. Probe tip electron delta-density. Used for electrostatic interaction when method is ‘fdbm’.

• A – float. Prefactor for Pauli repulsion when method is ‘fdbm’.

• B – float. Exponent used for Pauli repulsion when method is ‘fdbm’.

• fdbm_vdw_type – ‘D3’ or ‘LJ’. Type of vdW interaction to use when method is ‘fdbm’. ‘D3’ is for Grimme-D3 and ‘LJ’ uses standard Lennard-Jones vdW.

• d3_params – str or dict. Functional-specific scaling parameters for DFT-D3. Can be a str with the functional name or a dict with manually specified parameters. Used when method is ‘fdbm. See add_dftd3().

• lj_vdw_damp – int. Type of damping to use in vdw calculation when method is ‘fdbm’ and fdbm_vdw_type is ‘LJ’. -1: no damping, 0: constant, 1: R2, 2: R4, 3: invR4, 4: invR8.

• rot – np.ndarray of shape (3, 3). Rotation matrix applied to the atom coordinates.

• rot_center – np.ndarray of shape (3,). Point around which rotation is performed.

• local_size – tuple of a single int. Size of local work group on device.

• bRelease – Bool. Whether to delete data on device after computation is done.

• bCopy – Bool. Whether to copy the calculated forcefield field to host.

• bFinish – Bool. Whether to wait for execution to finish.

Returns:

FE – np.ndarray if bCopy==True or None otherwise. Calculated force field and energy.

prepareBuffers(atoms=None, cLJs=None, REAs=None, Zs=None, poss=None, bDirect=False, nz=20, pot=None, E_field=False, rho=None, rho_delta=None, rho_sample=None, minimize_memory=False)

Allocate all necessary buffers in device memory.

run_evalLJC_Hartree(FE=None, local_size=(32,), bCopy=True, bFinish=True)

Compute Lennard Jones force field at grid points and add to it the electrostatic force from an electric field precomputed from a Hartree potential.

Parameters:

• FE – np.ndarray or None. Array where output force field is copied to if bCopy == True. If None and bCopy == True, will be created automatically.

• local_size – tuple of a single int. Size of local work group on device.

• bCopy – Bool. Whether to return the calculated forcefield to host.

• bFinish – Bool. Whether to wait for execution to finish.

Returns:

FE – np.ndarray if bCopy == True or None otherwise. Calculated force field and energy.

run_evalLJC_QZs_noPos(FE=None, local_size=(32,), bCopy=True, bFinish=True)

Compute Lennard-Jones force field with several point-charges separated on the z-axis.

run_evalLJ_noPos(FE=None, local_size=(32,), bCopy=True, bFinish=True)

Compute Lennard-Jones forcefield without charges at grid points.

Parameters:

• FE – np.ndarray or None. Array where output force field is copied to if bCopy == True. If None and bCopy == True, will be created automatically.

• local_size – tuple of a single int. Size of local work group on device.

• bCopy – Bool. Whether to return the calculated electric field to host.

• bFinish – Bool. Whether to wait for execution to finish.

Returns:

FE – np.ndarray if bCopy == True or None otherwise. Calculated force field and energy.

run_gradPotentialGrid(pot=None, E_field=None, h=None, local_size=(32,), bCopy=True, bFinish=True)

Obtain electric field on the force field grid as the negative gradient of Hartree potential via centered difference.

Parameters:

• pot – HartreePotential or None. Hartree potential to differentiate. If None, has to be initialized beforehand with prepareBuffers.

• E_field – np.ndarray or None. Array where output electric field is copied to if bCopy == True. If None and bCopy == True, will be created automatically.

• h – float > 0.0 or None. Finite difference step size (one-sided) in angstroms. If None, the default value DEFAULT_FD_STEP is used.

• local_size – tuple of a single int. Size of local work group on device.

• bCopy – Bool. Whether to return the calculated electric field to host.

• bFinish – Bool. Whether to wait for execution to finish.

Returns:

E_field – np.ndarray if bCopy == True or None otherwise. Calculated electric field and potential.

setLvec(lvec, nDim=None)

setPP(Z_pp)

Set the atomic number of the probe particle. Required for calculating DFT-D3 parameters.

setQs(Qs=[100, -200, 100, 0], QZs=[0.1, 0, -0.1, 0])

tryReleaseBuffers()

Release all device buffers.

updateBuffers(atoms=None, cLJs=None, poss=None)

Update the content of device buffers.

class ppafm.ocl.field.HartreePotential(array, lvec, shape=None, ctx=None)

Bases: DataGrid

Sample Hartree potential. Units should be in Volts.

class ppafm.ocl.field.MultipoleTipDensity(lvec, shape, center=[0, 0, 0], sigma=0.71, multipole={'dz2': -0.1}, tilt=0.0, ctx=None)

Bases: TipDensity

Multipole probe tip charge density on a periodic grid.

Parameters:

• lvec – np.ndarray of shape (3, 3). Grid lattice vectors.

• shape – array-like of length 3. Grid shape.

• center – array-like of length 3. Center position of charge density in the grid.

• sigma – float. Width of charge distribution.

• multipole – Dict. Charge multipole types. The dict should contain float entries for at least of one the following ‘s’, ‘px’, ‘py’, ‘pz’, ‘dz2’, ‘dy2’, ‘dx2’, ‘dxy’ ‘dxz’, ‘dyz’. The tip charge density will be a linear combination of the specified multipole types with the specified weights.

• tilt – float. Tip charge tilt angle in radians.

• ctx – pyopencl.Context. OpenCL context for device buffer. Defaults to oclu.ctx.

class ppafm.ocl.field.TipDensity(array, lvec, shape=None, ctx=None)

Bases: DataGrid

Tip electron density. Units should be in e/Å^3.

interp_at(lvec_new, shape_new, array_out=None, local_size=(32,), queue=None)

Interpolate tip density onto a new grid. The tip is assumed to be cented on the origin, so the resizing of the grid happens in the middle of the grid and the corners remain fixed (the origin coordinates of the grids are ignored in the transformation).

Parameters:

• lvec_new – array-like of shape (4, 3). Unit cell boundaries for new grid.

• shape_new – array-like of length 3. New grid shape.

• array_out – pyopencl.Buffer or None. Output array. If None, then is created automatically.

• local_size – tuple of a single int. Size of local work group on device.

• queue – pyopencl.CommandQueue. OpenCL queue on which operation is performed. Defaults to oclu.queue.

Returns:

grid_out – TipDensity. New tip density grid.

subCores(xyzs, Zs, Rcore=0.7, valElDict=None)

Subtract core densities from the tip density.

Parameters:

• xyzs – np.ndarray of shape (n_atoms, 3). Coordinates of atoms.

• Zs – np.ndarray of shape (n_atoms,). Atomic numbers of atoms.

• Rcore – float. Width of core density distribution.

• valElDict – Dict or None. Dictionary of the number of valence electrons for elements. If None, then values in defaults.valelec_dict are used.

Returns:

TipDensity. New tip density with core densities subtracted.

ppafm.ocl.field.init(env)

ppafm.ocl.field.makeDivisibleUp(num, divisor)

ppafm.ocl.oclUtils

class ppafm.ocl.oclUtils.OCLEnvironment(i_platform=0)

Bases: object

loadProgram(fname)

printInfo()

printPlatformInfo()

updateBuffer(buff, cl_buff, access=<class 'pyopencl._cl.mem_flags'>)

ppafm.ocl.oclUtils.get_platforms()

ppafm.ocl.oclUtils.init_env(i_platform=0)

ppafm.ocl.oclUtils.print_platforms()

ppafm.ocl.relax

class ppafm.ocl.relax.RelaxedScanner

Bases: object

downloadPaths()

Get probe particle path array from device.

Returns:

paths – np.ndarray of shape scan_dim + (3,). xyz positions of probe particle at all scan points.

prepareAuxMapBuffers(bZMap=False, bFEmap=False, atoms=None)

prepareBuffers(FEin_np=None, lvec=None, FEin_cl=None, FEin_shape=None, scan_dim=None, nDimConv=None, nDimConvOut=None, bZMap=False, bFEmap=False, atoms=None)

prepareFEConv()

preparePosBasis(start=(-5.0, -5.0), end=(5.0, 5.0))

preparePossRot(pos0, avec, bvec)

releaseBuffers()

run(FEout=None, FEin=None, lvec=None, nz=None)

calculate force on relaxing probe particle approaching from top (z-direction)

run_convolveZ(FEconv=None, nz=None)

convolve 3D forcefield in FEout with 1D weight mask WZconv

run_getFEinStrokes(FEout=None, FEconv=None, FEin=None, lvec=None, nz=None, WZconv=None, bDoConv=False)

un-relaxed sampling of FE values from input Force-field (cl_ImgIn) store to cl_FEout

run_getFEinStrokesTilted(FEout=None, FEin=None, lvec=None, nz=None)

un-relaxed sampling of FE values from input Force-field (cl_ImgIn) store to cl_FEout operates in coordinates rotated by tipRot

run_getZisoFETilted(zMap=None, feMap=None, iso=0.0, nz=None)

get isosurface of input 3D field from given direction get map of 3D volume FE (e.g. electrostatic field) maped on 2D isosurface used to generate ElectrostaticMap returns zMap, feMap operates in coordinates rotated by tipRot

run_getZisoTilted(zMap=None, iso=0.0, nz=None)

get isosurface of input 3D field from given direction used to generate HeightMap operates in coordinates rotated by tipRot

run_izoZ(zMap=None, iso=0.0, nz=None)

get isosurface of input 3D field from top (z) used to generate HeightMap if cl_FEout is Forcefield it takes “z” where ( F(z) > iso )

run_relaxStrokesTilted(FEout=None, FEin=None, lvec=None, nz=None, bCopy=True, bFinish=True)

calculate force on relaxing probe particle approaching from particular direction

run_relaxStrokesTilted_convZ(FEconv=None, FEin=None, lvec=None, nz=None)

calculate force on relaxing probe particle approaching from particular direction

setScanRot(pos0, rot=None, zstep=None, tipR0=[0.0, 0.0, 4.0])

tryReleaseBuffers()

updateAtoms(atoms)

updateBuffers(FEin=None, lvec=None, WZconv=None)

updateFEin(FEin_cl, bFinish=False)

ppafm.ocl.relax.getInvCell(lvec)

ppafm.ocl.relax.init(env)

ppafm.ocl.relax.mat3x3to4f(M)

ppafm.ocl.relax.preparePoss(scan_dim, z0, start=(0.0, 0.0), end=(10.0, 10.0))

ppafm.ocl.relax.preparePossRot(scan_dim, pos0, avec, bvec, start=(-5.0, -5.0), end=(5.0, 5.0))

ppafm.ocl.relax.rotTip(rot, zstep, tipR0=[0.0, 0.0, 4.0])