Openmx manual

Please note, the manual was written by someone who was teaching themselves both R and OpenMx, and so some of the example scripts are cumbersome. A note on formatting I apologise for the erratic and inconsistent formatting in these blogposts.

Google's Blogger has many fine features, but control of formatting is not one of them. It has a mind of its own and will decide which font and spacing to use, despite instructions to the contrary. I don't want to spend hours delving in the html version to sort this out, so decided we'd have to live with it.

The scripts You may find it easier to copy the scripts from a Word document, rather than from this blog, in which case you can download the original document on which this blog is based from this site. Feedback Comments are enabled on this blog, and I'd welcome feedback and suggestions for improvement. I will try to incorporate any good ideas, but I can't provide technical support, as a I am not very expert myself, and b have a busy job that is only occasionally concerned with twin analysis.

I hope somebody out there finds this useful! Getting started in R. OpenMx is written in a programming language called R, which is free to download.

R is a powerful language for statistical computing, but much of the documentation is written for experts, and so it can be daunting for beginners. If you go to this website you will see instructions for how to download R.

Do not be put off by incomprehensible terms such as "CRAN mirror": just select a download site from the list provided that is geographically close to where you are.

You may then be offered further options that you may not fully understand. Just persevere by selecting the 'Windows' option from the 'Download and install R' section, and then select 'base', which at last takes you to a page with straightforward download instructions, and a link: ' Download R 2.

This will download an. I recommend accepting defaults. Installation of R will create a Start Menu item and an icon for R on your desktop.

Double clicking on the R icon starts the program. Starting R opens a window called R Console, in which you can type commands. This cursor will not be shown in the examples below, but it indicates that the console is awaiting input from you.

This starts a web-based interface to on-line help. You may want to briefly explore this before going further. As with other programming languages, you hit Enter at the end of each command.

Just to familiarise yourself with the console, type:. R evaluates the expression and you see output:. The [1] at the beginning of the output line indicates that the answer is the first element of a vector. All these four values are given in units eV relative to Fermi level. The inner window should be fully inside of the outer window.

If the two boundaries of inner window are equal to each other, it means inner window is not defined and not used in calculation. There is no default values for outer window, while 0.

Bottom Top 0. Bottom 0. If you want to restart the minimization of MLWFs by reading the overlap matrix elements from files, the outer window should not be larger than that used for calculating the stored overlap matrix. Either equal or smaller is allowed. The inner window can be varied within the outer window as you like when the restart calculation is performed.

This would benefit the restarting calculation or checking the dependence of MLWFs on the size of both the windows. If the initial guess is required, a set of local functions with the same number of target MLWFs should be defined. Bloch wave functions inside the outer window will be projected on to them.

Therefore, these local functions are also called as projectors. The following steps are required to specify a projector. Define local functions for projectors Since the pseudo-atomic orbitals are used for projectors, the specification of them is the same as for the basis functions.

An example setting, for silicon in diamond structure, is as following: Species. Species Si Si7. In fact, the pseudopotential de- fined in this line will not be used. It is given just for keeping the consistence of inputting data structure. One can use any PAO as projector. Also the use of only a single basis set is allowed for each l-component.

Projectors proj1-sp3 0. The next two sets of three numbers define the z-axis and x-axis of the local coordinate system, respectively, where each axis is specified by the vector defined by three components in xyz-coordinate. In the second line the local axes are the same as the original coordinate system. The orbital used as projector can be the original PAOs or any hybrid of them.

A list of supported PAOs and hybridizations among them can be found in Table 6. Any name other than those listed is not allowed. The projector can be centered anywhere inside the unit cell.

There is no default setting for it. To use finite difference approach for calculating k-space differentials, b vectors connecting neighboring k points are searched shell by shell according to the distance from a central k point.

Default value is 12 and it should be increased if failure in finding a set of proper b vectors. The problem may happen in case of a system having a large aspect ratio among unit vectors, and in this case you will see an error message, while the value 12 works well in most cases. The hybridization is done within the new coordinate system defined by z-axis and x-axis. MaxShells 12 default value is The first step is to minimize the gauge invariant part of spread function by disentangling the non-isolated bands.

The second step is the same as isolated band case [73]. The gauge dependent part is optimized by unitary transformation of the selected Bloch wave functions according to the gradient of spread function. For the first step, three parameters are used to control the self-consistence loop. They are the maximum number of SCF loops, the convergence criterion, and the parameter to control the mixing of input and output subspace projectors, respectively.

Steps default Wannier. Criterion 1e default 1e-8 Wannier. Para 0. One is a steepest decent SD method, and the second one is a conjugate gradient CG method.

The third one is a hybrid method which uses the SD method firstly and then switches to the CG method. In the CG method, a secant method is used to determine the optimized step length. StepLength 2. Steps 5 default 5 Wannier. The code will read the overlap matrix as well as the eigenenergies and states from the disk file.

One should keep in mind that the outer window and k grid should be the same as those used for calculating the stored overlap matrix and eigenvalues. Consistence will be checked in the code.

The inner window, initial guess of MLWF as well as the convergence criteria can be adjusted for restarting optimization.

As an example, the interpolated band structure of Si in diamond structure is shown together with its original band structure in Fig. Plot on default off Wannier. Users can set the supercell size for plotting MLWF. Figure 34 b shows one of the eight converged MLWFs from four valence states and four conduction states near Fermi level of Si in diamond structure. It is obtained with an initial guess of sp3 hybrid.

To monitor the optimiza- tion progress, the following method may be helpful. The following example is for Si. SD 4. CG 58 1. The details of the equation can be found in Ref. SD CG 57 The initial guess is sp3 hybrids. The initial guess is sp2 hybrids and pz orbitals on carbon atoms. They have different extension k,b names. More detailed information of the four files will be given below.

The first line of this file is the description of the numbers in the second line. The numbers from left to right in the second line are the number Nwin of included bands within the outer window, the number of k points, the number of b vectors, k,b the number of spin component, respectively.

The next lines are data blocks of Mmn. The most outer loop is for spin component. The next is the loop of k points and then b vectors. The most inner loops are the band index n and m, respectively. In each block, the first line are 5 numbers. From the second line are the real and imaginary part of each matrix element.

The first line of the file is the description of the whole file. Obviously, the four numbers in the second line are the number Nwin of bands within the outer window, the number of k points, the number of target MLWFs and the number of spin component, respectively.

Similarly, the data blocks are written in loops. The most outer loop is spin component and then k points, target MLWFs and number of bands. As described in the first line of this file. In each block, the first three integers are the band index, the index of MLWFs and index of k points, respectively. The next are real and imaginary of that matrix element. Next is m n k and elements. Spin is the most outer loop. The first line is the Fermi level of system.

The number of bands is indicated in the second line of the file. The next data are mainly in two parts. The first part is the eigenenergies and the second one is the corresponding eigenstates. In each part, the loop of spin component is the most outer one.

The next loop is k points, followed by band index. For eigenstates, there is one more inner loop for the basis set. WF kpt 1 0. The number of MLWFs, number of lattice vectors inside of Wigner-Seitz supercell are in the second and third line, respectively.

The unit cell vectors are given in the fifth, sixth and seventh lines. Spin polarization, whether it is a non-spin polarized calculation or a spin polarized one with collinear or noncollinear magnetic configuration, is given in the eighth line. The ninth line gives the Fermi level. From the tenth line, the block data starts. The outer most loop is spin component. The next loop is for R and the last two are loops of m and n, respectively. Each R is written at the first line of each block together with its degeneracy.

The index of m and n is printed and followed by the real and imaginary parts of hopping integrals in each line. The reference results were calculated using a Xeon cluster machine. Unlike O N methods developed so far the approach is a numerically exact alternative to conventional O N 3 diagonalization schemes in spite of the low-order scaling, and can be applicable to not only insulating but also metallic systems in a single framework.

The well separated data structure is suitable for the massively parallel computation as shown in Fig. However, the advantage of the method can be obtained only when a large number of cores are used for parallelization, since the prefactor of computational efforts can be large.

When you calculate low-dimensional large-scale systems using a large number of cores, the method can be a proper choice. The electric temperature of K and 80 poles for the contour integration are used. For comparison, the speed-up ratio for the parallel computation of the conventional scheme using Householder and QR methods is also shown for the case with a single thread. The elapsed time at cases pointed by arrow is also shown for both the low-order scaling and conventional methods.

The input file is C60 LO. Method Total energy Hartree Computational time sec. Low-order Thus, it is possible to perform geometry optimization. However, calculations of density of states and wave functions are not supported yet. The number of poles in the contour integration [51] is controled by a keyword: scf. ON2 90 The number of poles to achieve convergence does not depend on the size of system [75], but depends on the spectrum radius of system.

If the electronic temperature more K is used, poles is enough to get sufficient convergence for the total energy and forces. We expect that the crossing point between the low-order scaling and the conventional methods with respect to computational time is located at around atoms when using more than cores for the parallel computation, although it depends on the dimensionality of system.

In this method, a 2-dimentional periodic and 1- dimentional optional boundary conditions are imposed on a model cell Fig. An isolated slab, charged slab, and a slab under an uniform electric field can be treated by introducing the following combinations of semi-infinite media ESMs.

An isolated slab model can be used for investigations of a polarized substrate, and charged slab model is applicable to a simulation of an electrode surface. A slab model under an electric filed sandwiched between two ideal-metal media would be appropriate for a material located in a metal capacitor. The slab and ESMs are placed parallel to the y-z plane.

The a-axis of the cell is perpendicular to the b-c plane and is parallel to the x-axis. Two periodic boundary conditions are set in y- and z-axis directions 3.

The origin of the x-axis is set at the cell boundary. A fractional coordinate for x-axis is designated between 0 and 1. In this case, note that the total charge of a calculation system should be neutral. One can deal with charged systems.

By using the following keyword, one can impose a uniform electric field on a calculation system; ESM. The electric filed is decided by the length of the cell, a, and the potential difference. A surface-model slab and a liquid region should be located as shown in Fig. In order to restrict liquid molecules within a given region, an cubic barrier potential can be introduced by using the following keyword see Fig.

The number of doped charge is Each plot is obtained as a difference in difference charge or difference Hartree potential with reference to a neutral slab with the same ESMs. It is also recommended to fix positions of atoms on the bottom of a surface-model slab during a MD run.

It can be seen that segregation of the doped charge in the slab happened due to the attractive interaction between the doped and the corresponding mirror charges. Geometry optimization of a precursor 2. Geometry optimization of a product 3. Optimization of a minimum energy path MEP connecting the precursor and product where in the three calculations users have to keep the same computational parameters such as unit cell, cutoff energy, basis functions, pseudopotentials, and electronic temperatures to avoid numerical inconsistency.

The spring constant is given by MD. Const 0. Some of them are listed below: c2h4. The second column is the total energy of each image. The third and fourth columns are interval bohr between two neiboring images and the distance borh from the precursor in geometrical phase space.

Name and is a serial number for each image, is also generated, since each calculation for each image is basically done as an indenpendent OpenMX calculation with a different input file. After the successful calculation, you may get the history of optimization and change of total energy along MEP as shown in Fig. In such a case, one has to continue the optimization as a new job starting from the last optimization step in the previous job.

The file contains a series of atomic coordinates for images in the last step. However, in some cases the geometrical structure of images generated on the straight line can be very erratic so that distance between atoms can be too close to each other.

In this case, one should explicitly provide the atomic coordinates of images. The user defined initial path can be provided by the same way as for the restarting. SpeciesAndCoordinates 1 Si Images, the atomic coordinates need to be provided. Also, it is required for a keyword to be switched on as scf.

However, there is no guarantee that the SCF iteration con- verges for all the images. In order to monitor the SCF convergence for all the images, temporary files can be checked by users. However, it would be better to provide a proper number of processes for the MPI parallelization which can be divisible by the number of images given by MD.

Images, in order to achieve a good load balance in the MPI parallelization. It is noted that the number of processes for the MPI parallelization can exceed the number of atoms unlike the conventional calculation. The translation tends to generate a confusing movie in the visualization of the result. Only three routines are added to implement the NEB functionality. They are neb.

The main routine is neb. It may be easy to implement related methods in neb. The method is nothing but calculation of partial charge density in an energy window measured from the chemical level.

The calculation of the partial charge density is performed by the following keywords: partial. Since the calculation of the partial charge density is performed during calculation of the density of states DOS , the following keywords have to be specified as well: Dos. As an example, a simulated STM image of a graphene layer is shown in Fig. Figure Simulated STM image of a graphene layer, where partial.

The following keywords are relevant to the DFT-D2 method. Unit DFTD. The scaling factor in Eq. By considering the periodicity or non-periodicity of each atom, the interaction is automatically cut when they are non-periodic.

Since OpenMX uses localized orbitals as basis function, it is very important to take account of basis set superposition error BSSE when we investigate an effect of a weak interaction such as vdW interaction.

Step 0. Figure Total energy vs. When you perform calculations of the density of states by the following keywords: Dos.

Name, and the file can be visualized by XCrySDen [56]. It is noted that a large number of k-points should be used in order to obtain a smooth Fermi surface. As an example, Fermi surfaces of the fcc Ca bulk are shown in Fig. Since two sorts of bands intersect with the Fermi energy chemical potential , two Fermi surfaces are shown in a and b.

If you analyze the difference between two states, this tool would be useful. Thus, you can easily visualize the difference using many softwares, such XCrySDen [56] and Molekel [55]. In fact, Fig. If you analyze the difference between two geometries, this tool would be useful. These figures were visualized by XCrySDen. The calculated result appears in the standard output your display. Bond lengths below the cutoff bond length are taken into account for the RMSD calculation. Figure 43 shows vectors corresponding to the deviation of atomic coordinates in optimized structures and the difference of total charge density between a neutral and one electron doped glycine molecule.

We see that the large structural change seems to take place together with the large charge deviation. This example illustrates that the tool would be useful when we want to know how the structure is changed by the charge doping and the electric field. Also, the coordinates of the system A must be the same as in the calculation i.

To use the same origin as in the calculation i rather than the use of an automatically determined origin, you have to include the following keyword in your input file: scf. Then, you will have a cube file for charge spin density. Also, the coordinates of the system B must be the same as in the calculation i.

Therefore, the use of the minimum cell size is desirable in terms of computational efficiency. OpenMX supports the requirement. However, we see a platform dependency of lapack routines to solve the tridiagonalized matrix with respect to computational robustness.

So, three different lapack routines are available in OpenMX Ver. For further details, see the lapack website [85]. If you want to use the Kohn-Sham Hamiltonian, the overlap, and the density matrices, Then these data can be utilized by the following steps. Name in your input file. If any program bug is introduced, they will not be consistent with each other.

However, the dynamic memory allocation causes often a serious memory leak which wastes the memory used as the MD steps increase. In such a case, you may observe erratic eigenvalues. To avoid the overcompleteness, a small number of optimized basis functions should be used. In such a case, one has to mix the charge density very slowly, indicating that the number of SCF steps to get the convergence becomes large unfortunately. It is expected that the forum is utilized for sharing tips in use of OpenMX and for further code development.

Moreover, the author, Taisuke Ozaki, possesses the copyright of the original version of this program package. We cannot offer any guarantee in your use of this program package.

However, when you report program bugs, we will cooperate and work well as much as possible together with you to remove the problems. Acknowledgment One of us T. One of us T.

References [1] P. Hohenberg and W. Kohn, Phys. Kohn and L. Sham, Phys. Ceperley and B. Alder, Phys. Perdew and A. Zunger, Phys. B 23, Perdew and Y. Wang, Phys. B 45, Perdew, K. Burke, and M. Ernzerhof, Phys. Barth and L. Hedin, J. C: Solid State Phys. Sticht, and A. Williams, J. F: Met. Sticht, K-H. Matter 1, Oda, A. Pasquarello, and R. Car, Phys. MacDonald and S. Vosko, J. Kurz, F. Forster, L. Nordstrom, G, Bihlmayer, and S. Blugel, Phys. B 69, King-Smith and D. This is typical example of matrix keyword.

User can specify explicitly this argument using python list object. See the official OpenMX manual for more detail. However, single node calculating is not a good way to run heavy DFT calculation. Parallel computation is inevitable. In OpenMX calculator, user may choose the way to run. There are two ways to excute the code. MPI method can be applied in general. To use it, put the mpi dictionary as a kwargs. It will define the name of file.

Restart the calculation. Dictionary contains the key processes and threads. You should probably consider just deleting unwanted variables from the file and then resaving in. An alternative is to save in. The downside of this is that blank entries will not be treated as missing, and so you need to do some extra processing to avoid a file read error. There are probably better ways of doing this, but this is what worked for me!

Open the. How you do this depends on how missing data are represented in your file. You may have used code numbers such as , or the missing data may be represented by blanks.