Merge pull request #3 from nanotech-empa/release/1.4.0

Release/1.4.0

.gitignore

@@ -1,3 +1,4 @@
 .ipynb_checkpoints/
 __pycache__/
+.vscode
 

README.md

@@ -11,6 +11,17 @@ Features include:
 * Orbital hybridization analysis for adsorbed systems
 * High-resolution STM (HRSTM) simulations
 
+Requirements:
+* numerical and scientific python libraries `numpy`, `scipy`
+* Atomistic simulation environment `ase`
+* `mpi4py` provides MPI parallelization
+
+When everything is set up correctly, the bash scripts in `examples/` folder can be executed without any further input and illustrate the usage of the various scripts. For example `example/benzene_stm/run_stm_sts_from_wfn.sh` evaluates the STM/STS signatures of isolated benzene at each orbital energy (`out/orb/`) as well as in an arbitrary energy range (`out/stm/`). The corresponding CP2K calculation is included in the repository.
+
+**NB: In all cases, the underlying DFT calculation has to be performed with the diagonalization algorithm rather than orbital transformation (OT).**
+
+### Evaluating molecular orbitals on an arbitrary grid
+
 Most of the functionality of this library is built on top of the possibility to evaluate the Kohn-Sham orbitals encoded in the `.wfn` file on an arbitrarily defined grid. This is illustrated by the following script applied for a nanographene adsorbed on a Au(111) slab (total of 1252 atoms and 10512 electrons):
 
 ```python
@@ -25,16 +36,16 @@ cgo.load_restart_wfn_file("./PROJ-RESTART.wfn", n_occ=2, n_virt=2)
 
 ### Evaluate the orbitals in the specific region ###
 cgo.calc_morbs_in_region(
-    dr_guess = 0.15,                        # grid spacing
-    x_eval_region = None,                   # take whole cell in x
-    y_eval_region = [0.0, cgo.cell[1]/2],   # half cell in y
-    z_eval_region = [36.0, 44.0],           # around the molecule in z
+    dr_guess = 0.15,                          # grid spacing
+    x_eval_region = None,                     # take whole cell in x
+    y_eval_region = [0.0, cgo.cell_ang[1]/2], # half cell in y
+    z_eval_region = [19.0, 24.0],             # around the molecule in z
 )
 
 cgo.write_cube("./homo.cube", orbital_nr=0)
 ```
 
 The evaluated HOMO orbital in the defined region:
 
-
 <img src="examples/example.png" width="600">
+

cp2k_spm_tools/__init__.py

@@ -1,3 +1,5 @@
-"""Atomistic simulation tools (AST)
+"""CP2K Scanning Probe Microscopy tools
 
 """
+
+__version__ = 1.4.0

cp2k_spm_tools/common.py

@@ -57,7 +57,6 @@ def parse_eval_region_input(eval_reg_inp, ase_atoms, cell):
                     print("Unrecognized option ", eval_regions_inp[i][j])
                     sys.exit(1)
                 eval_regions[i][j] = float(reg_str)
-                eval_regions[i][j] *= ang_2_bohr
 
             else:
                 ref_at_pos = reg_str[0]
@@ -74,7 +73,6 @@ def parse_eval_region_input(eval_reg_inp, ase_atoms, cell):
                 eval_regions[i][j] = (np.min(sel_positions[:, i]) + ref_shift_val
                                         if ref_at_pos == 'n' else
                                     np.max(sel_positions[:, i]) + ref_shift_val)
-                eval_regions[i][j] *= ang_2_bohr
 
         if np.abs(eval_regions[i][0] - eval_regions[i][1]) < 1e-3:
             eval_regions[i][0] = eval_regions[i][1]

cp2k_spm_tools/cp2k_grid_orbitals.py

@@ -24,7 +24,7 @@
 
 from mpi4py import MPI
 
-ang_2_bohr = 1.0/0.52917721067
+ang_2_bohr = 1.88972612463
 hart_2_ev = 27.21138602
 
 class Cp2kGridOrbitals:
@@ -72,12 +72,15 @@ def __init__(self, mpi_rank=0, mpi_size=1, mpi_comm=None, single_precision=True)
         # Orbitals on discrete grid
         self.morb_grids = None
         self.dv = None # [dx, dy, dz] in [au]
-        self.origin = None
+        self.origin = None # origin point of the evaluation grid
         self.eval_cell = None
         self.eval_cell_n = None
 
         self.last_calc_iz = None # last directly calculated z plane (others extrapolated)
 
+    @property
+    def cell_ang(self):
+        return self.cell/ang_2_bohr
 
     ### -----------------------------------------
     ### General cp2k routines
@@ -297,7 +300,7 @@ def read_basis_functions(self, basis_set_file):
 
     def load_restart_wfn_file(self, restart_file, emin=None, emax=None, n_occ=None, n_virt=None):
         """
-        Reads the specified molecular orbitals from cp2k restart wavefunction file
+        Reads the specified molecular orbitals from cp2k restart wavefunction file.
         If both, energy limits and counts are given, then the extreme is used
         Note that the energy range is in eV and with respect to HOMO energy.
         """
@@ -519,7 +522,7 @@ def calc_morbs_in_region(self, dr_guess,
         Puts the molecular orbitals onto a specified grid
         Arguments:
         dr_guess -- spatial discretization step [ang], real value will change for every axis due to rounding  
-        x_eval_region -- x evaluation (min, max) in [au]. If min == max, then evaluation only works on a plane.
+        x_eval_region -- x evaluation (min, max) in [ang]. If min == max, then evaluation only works on a plane.
                         If left at None, the whole range of the cell is taken.
         pbc -- determines if periodic boundary conditions are applied in direction (x, y, z) (based on global cell)
         eval_cutoff -- cutoff in [ang] for orbital evaluation
@@ -531,12 +534,14 @@ def calc_morbs_in_region(self, dr_guess,
         eval_cutoff *= ang_2_bohr
         reserve_extrap *= ang_2_bohr
 
+        eval_regions_ang = [x_eval_region, y_eval_region, z_eval_region]
+        eval_regions = [np.array(er)*ang_2_bohr if er is not None else er for er in eval_regions_ang]
+
         global_cell_n = (np.round(self.cell/dr_guess)).astype(int)
         self.dv = self.cell / global_cell_n
 
         ### ----------------------------------------
         ### Define evaluation grid
-        eval_regions = [x_eval_region, y_eval_region, z_eval_region]
         self.eval_cell_n = np.zeros(3, dtype=int)
         self.origin = np.zeros(3)
 
@@ -764,6 +769,8 @@ def extrapolate_morbs_spin(self, ispin, vacuum_pot=None, hart_plane=None, use_we
 
         print("Extrapolation time: %.3f s"%(time.time()-time1))
 
+
+
     ### -----------------------------------------
     ### Export data
     ### -----------------------------------------
@@ -775,7 +782,7 @@ def write_cube(self, filename, orbital_nr, spin=0, square=False):
             print("R%d/%d is writing HOMO%+d cube" %(self.mpi_rank, self.mpi_size, orbital_nr))
 
             energy = self.morb_energies[spin][local_ind]
-            comment = "E=%.8f eV (wrt HOMO)" % energy
+            comment = "E=%.8f eV (wrt middle of gap)" % energy
 
             if not square:
                 c = Cube(title="HOMO%+d"%orbital_nr, comment=comment, ase_atoms=self.ase_atoms,