MAPXCH : Import theoretical electron density maps

Authors: Doug du Boulay, Materials and Structures Laboratory, Tokyo Institute of Technology, Nagatsuta, Midori-Ku, Yokohama Japan.

MAPXCH currently imports a abinitio (DFT) 3D electron density map from the open source ABINITproject or from Wien2klapw5-3D into an Xtal archive file.

Purpose

To compare theoretical and experimental charge densities it is useful to know that all densities are on the same scale, with the same contour intervals and the same map orientations, so that topographic features can be reliably compared. To this end MAPXCH was written to import the electron density maps of the completely unrelated and unaffiliated external density functional theory programs Wien2k and ABINIT into Xtal, to exploit the relatively reliable contouring machinary that already exists herein (amongst other things).

Using Wien2k data

At last checking, the Wien2kDFT code included one program, lapw5, intended for determining the electron density on a planar grid of values for contouring in external software. Florent Boucherkindly made available (to existing Wien2k users) a modified version of lapw5, heretofore named lapw53D, with which it was possible to compute the electron density on a 3D grid of values, in much the same manner as an existng planar lapw5 calculation.

The program lapw53D requires a default filename mapping file - lapw53D.defwith the contents:

5 ,'case.in53D', 'old',    'formatted',0
6 ,'case.output5',   'unknown','formatted',0
8 ,'case.struct',    'old',    'formatted',0
9 ,'case.clmval',    'old',    'formatted',0
10,'case.3dmap',     'unknown','unformatted',0
11,'case.clmvaldn',  'unknown','formatted',0
12,'case.sigma',     'unknown','formatted',0
20,'case.rho_onedim','unknown','formatted',0
21,'case.rho',       'unknown','formatted',0

For use with Xtal we need a lapw53D control file case.in53Dwith the contents:

           0           0           0      10000000
    10000000           0           0      10000000
           0    10000000           0      10000000
           0           0    10000000      10000000
3 2 3
97 97 49
DIFFADD
ANG VAL NODEBUG
NONORTHO

Note especially the grid points 97 97 49 which will vary on a study by study basis, but if they are too large, the whole map cannot be squeezed into the F77 fixed length Xtal memory. For Xtal it is necessary that the unduplicated grid dimensions be divisible by 2, 3 and/or 5 and for trigonal symmetrys divisibility by 3 is essential.

So, having run your Wien2k DFT SCF calculations to convergence, run lapw53D to calculate the 3D grid data points which are stored in the raw binary case.3dmapfile.

Then, run Xtal with the command MAPXCH w2k2map inpfile case to read both the case.structfile to obtain unit cell, symmetry and atom information and, also the case.3dmapto transform the 3D gridded Wien2k density map into the Xtal map file as a full cell map. After the Xtal archive is created STARTX and ADDATM are run in update mode to check the symmetry and atom site information, then the archive is coppied to the .mapfile.

If you don't actually have a case.3dmapfile, the MAPXCH w2k2map command should still read in the case.struct, which should be sufficient to, for instance, run PIG and produce an equivalent .pdbfile, or to run CIFIO cifout to generate a CIF.

Using ABINIT data

To obtain a valence density output from the program ABINIT it is sufficient to add the command

prtden 1

to your ABINIT control file. Then at succesful completion of the job all cell, symmetry and atom site information is written, along with the resultant valence density, to the file case_DEN. In contrast to Xtal, ABINIT plays hard and fast with the symmetry and cell in order to minimize cell volumes and optimize computation times. As a result some deductive logic is required to transform the 3D map data and cell to a more standard crystallographic setting. This can be lossy and the MAPXCH code is experimental.

One problem with the ABINIT case_DENfile is that the entire electron density map is written as a single Fortran90 record of double precision data. Reading it within Xtal stresses the memory limitations, severely limiting the size of map that can be accomodated.

As with the Wien2k data input method, after creating the Xtal archive it is passed to STARTX and ADDATM to check the symmetry and atom sites, before being coppied to a .mapfile for contouring within xtal.

Using VALRAY data

It turned out to be impossible to read the 3D electron density from VALRAY, because it doesnt create such a file. Even the 2D data optionally stored for a slice in the bindaffile comes with no parameters describing its origin, orientation, dimensions etc., which is a shame. Nevertheless you can read in all the reflection and phase information from the bindaffile and use it to calculate a residual Fourier density map, in Xtal. No multipoles, but thats life.

Using XD data

Having run an XD multipole refinement with xdlsm, it is then necessary to run xdfour to evaluate the electron density on an appropriate three dimensional grid of data points as shown in the examples to follow. MAPXCH can then read in the unit cell and symmetry from the xd.masfile, the refined coordinates from the xd.resfile and the 3d grid from xd_fou.grd.

The only real reason for bothering with this is that the atom site labeling in CONTRS is essentially automatic, the postscript rendering from PREVUE con or PLOTX con works ok under linux and, thereafter SURFIN gives a relatively quick and nasty interactive view of the XD electron density.

File Assignments

  • Optionally read ABINIT case_DENfile

  • Optionally read Wien2k case.structand case.3dmapfiles

  • Initializes an xtal archive bdf and creates an xtal mapfile.

Examples

MAPXCH w2k2map inpfile P63 :  read Wien2k P63.tmp map file
UOV 0.22

Read a Wien2k 3D map file, presumably a difference electron density map, straight into an xtal .map file

MAPXCH abi2map noterm abifile sin_o_DEN :  read ABINIT map
UOV 0.22                   : spawns STARTX and ADDATM
ADDATM upd:                : reoveride builtin value
UOV 0.00

ADDREF ffac :              generate dummy reflection data
reduce nocon
limits *9 yes :            whole sphere
hklgen hkl  frel sigf

SORTRF aver 1 frel :       remove redundant reflections

RFOURR : pr 100 :          reverse fourier transform valence
density

MODHKL :                   purge unphased reflections
purge -4.+19 $1 1 1801

FC frt noex nod  list  :   Determine Frel as
Sqrt((Acore+Aval)^2+(Bcore+Bval)^2)

CRYLSQ ov cy 10 wu :       optimize overall temp
noref skf
noref Si
noref N

FOURR  :                   Calculate a difference map
grid 60 60 24

This example reads an ABINIT valence electron density map into an xtal input .aanarchive file for updating using STARTX and ADDATM before copying to the standard .mapfile. Dummy reflections are determined (assuming λ=0.5Ĺ) and a reverse Fourier transform calculates the phase as well as A and B. The Fc frt option calculates an effective Frel by adding a pseudo core to the valence electrons [291]. Then total density, promolecular structure factors are computed in CRYLSQ and a difference Fourier map is evaluated.

MAPXCH abi2map noterm abifile sin_o_DEN :  read ABINIT map
UOV 0.00

ADDATM upd     
UOV 0.000

STARTX upd     :                Add valence density form factors
celcon Si
celcon N
setid formfx
Si 0.  4.
Si 0.00836484041  3.99153447
: ... truncated for brevity. Valence form factors from SCATOM
N  1.98792839  4.64325094
setid
exper 1  0.9 :                  change the wavelength
end

ADDREF ffac :                   generate dummy reflection data
reduce nocon
limits *4 1.0   *9 yes :        whole sphere
hklgen hkl  frel sigf

SORTRF aver 1 frel :            remove redundant reflections

RFOURR th 0.000001 : pr 100 :   reverse fourier transform valence
density
idnums $1 304 801 802 700   :

MODHKL :                        purge unphased reflections
purge -4.+19 $1 1 1801
FC noex nod list :              calculate Fc for pure valence
electrons

FOURR  :                        Calculate a valence density
difference map
grid 60 60 24

Calculate a true ABINIT difference valence density map.

A VALRAY import example

COMPID ximpor
MAPXCH val2map inpfile bindaf.in
FOURR pset 1 full             : Use VALRAY Fcal and phase     
grid 144 60 48
TITLE  Figure 1.
SLANT file nprint:            : compute a residual density slice
...
...

An XD import example, purely for electron density contour maps.

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   MODULE *XDFOUR
SELECT   fobs *fmod1 *fmod2 print snlmin 0. snlmax 2.
GRID  3-points perp *cryst
LIMITS xmin 0.0 xmax 1.0 nx 144
LIMITS ymin 0.0 ymax 1.0 ny 60
LIMITS zmin 0.0 zmax 1.0 nz 48
   END XDFOUR
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

COMPID ximpor
:   read in  xd.mas, xd.res and xd_fou.grd
MAPXCH xd2map inpfile xd
:                             
SLANT file nprint:            : compute a fmod1-fmod2 slice ...
...


[291] Don't do this. It assumes the missing core electron density is something resembling an inert gas atom which is a very bad approximation.