EPIC is a front end to SCEPIC, which helps to generate inputs, and optionally execute SCEPIC.
Prerequisite software: Molcas (8.6) or OpenMolcas.
Please, verify that Molcas or OpenMolcas is installed and fully functional (by running the tests). Use the corresponding documentation to make one of this installations. For Molcas installation the driver 'molcas' should be in the PATH, for OpenMolcas 'pymolcas' should be in the PATH.
SCEPIC is a collection of Julia functions (available at https://gitlab.com/ErnDLar/scepic or molcas.org/SCEPIC) that requires no real installation.
Under Linux Julia interpreter can be installed bu standard package manager, e.g. for Ubuntu by sudo apt install julia. Julia interpreter can be also installed from the source (https://julialang.org/downloads/), along with the (standard) libraries: Base.Threads, LinearAlgebra, Printf. So far, the SCEPIC has only been tested on macOS (Big Sur) and Linux (Ubuntu, CentOS).
In order to run SCEPIC, some environment variables has to be set:
export SCEPIC_ROOT=$PWD #it is for example, /home/valera/scepic. export JULIA_LOAD_PATH=$SCEPIC_ROOT/src/embedding:$SCEPIC_ROOT/src/scepic:$JULIA_LOAD_PATH export MOLCAS=/home/valera/OpenMolcas/build/ #location of Molcas or OpenMolcas export SCEPIC_SCRATCH=/tmp/scepic #SCEPIC_SCRATCH directory will be used by Molcas/OpenMolcas for # the storage of temporary files (WorkDir). # Note that the directory should exist.
Create a directory and copy files from $SCEPIC_ROOT/tests/embedding/in/ to it. run julia driver.jl >log. (it can take some time). The resulting files can be compared with files, located at $SCEPIC_ROOT/tests/embedding/out/. project.xfield file can be used directly in Molcas calculation, project.inp is a prototype of an input file for Molcas and project_xfield.xyz file can be used for visualization of point charges. Note that the atoms in project_xfield.xyz are Na for cations and F for anions.
Create a directory and copy files from $SCEPIC_ROOT/tests/scepic/in to it. run julia driver.jl >log. (it can take some time). The resulting files can be compared with files, located at $SCEPIC_ROOT/tests/scepic/out/
Note that computed files might be slightly different depending on the version of Molcas you use.
In some trivial cases it is enough to use EPIC script, which generates all necessary files with default options. In this section we will describe how this script works and later the more detailed explanation of the parts of SCEPIC will be given.
You need to create an XYZ file for your periodic structure, with an extended format. The file should contain all atoms in the unit cell, but in addition to that,
#Origin 0 0 0 #Vector1 2.0 0 0 #Vector2 0 2.0 0 #Vector3 0 0 2.0
Mg 0 0 0 +2 sp
In case if one atom located at different Wyckoff positions, the element name must be labelled, by adding underscore sign and any label, e.g.
Mg_a 0 0 0 +2 sp
After creating the file, execute ./EPIC file.xyz. The script will create file runme, which should be executed and at the end a file, named ECP is created.
Before running runme file you may alter the driver scripts.
Note, that the version of EPIC from 231127 not always correctly handles some non cubic cells. Be sure that XYZ coordinates can be converted into partial coordinates in a way that the partial coordinates are in [0,1] interval. So, in teh example above, a position -2.0 0 0 might be problematic.
EPIC understands the flags from the command line, and/or from the comment line in the XYZ file (at the later case, the flags must be preseeded by -EPIC). A command EPIC file.xyz -aimp=20 is the same as command EPIC file.xyz if the second line in file.xyz file looks like
comment line -EPIC -aimp=20
The most important flag is -cluster=cluster.xyz file. It tells EPIC to generate a prototype of input for Molcas, instead of creating inputs for AIMP generation.
Another flag is -run which tells EPIC to run SCEPIC. Without this flag user must run SCEPIC via generated files (runme, or runmec (for the cluster mode)).
The full list of EPIC flags can be found it its source code.
The program works in units of Ångstrom (Å). Unit-cell information has to be specified in an atomic simulation environment-extended xyz-file. In principle, no more information than the unit-cell is necessary to create a nanocluster to represent an ionic crystal with the GPEE method. To make the program more user-friendly, an additional file specifying the desired QM-region can be given. A Molcas-type input will then be prepared, where the ions included in the QM-region have been excluded from the formal point-charges generated by the embedding program.
Since the method is intended to function with AIMPs, keywords associated with automatic replacement of point-charges with AIMPs are also present. The only effect from these keywords it to automatically write the Molcas input-file with AIMP specification present. Currently, there is no method in the code to distinguish direct-border AIMPs from the remaining AIMPs. Therefore, if the addition of extra basis function on direct-border AIMPs is desired to improve cluster-host orthogonality, these should be included in the xyz-file for the QM-region.
An example of how an input script to embedding program might look is given below. This example is also present in the test-directory.
See all allowed keywords in Table ![[*]](crossref.png)
.
Remember to add the embedding source directory to the Julia load path before running it as "export JULIA_LOAD_PATH=path_to_embedding_source:$JULIA_LOAD_PATH".
After setting the load path, the script is execute on the command-line as "julia name_of_script.jl".
Example of an atomic simulation environment-extended xyz-file, specifying a MgO-unit cell. For the purpose of the embedding program, only the Lattice="..." part of the comment-line is important.
8 Lattice="4.256484 0.0 0.0 0.0 4.256484 0.0 0.0 0.0 4.256484" Mg 0.00000000 0.00000000 0.00000000 ... O 2.12824200 0.00000000 0.00000000 ...
Example of a simple script to generate a nanocluster for MgO. This example is also present in the test-directory.
using embedding_factory
uc_file = "./MgO_bulk.xyz"
qm_file = "./MgO_clust.xyz"
charge_dict = Dict("Mg" => 2, "O" => -2)
input = Dict("elmoment" => 4,
"unit_cell_file" => uc_file,
"qm_file" => qm_file,
"cluster_radius" => 40,
"aimp_size" => 3,
"charge_dict" => charge_dict)
embedding(input)
A quick summary of SCEI: The ions of a material, e.g., MgO are embedded in a host of AIMPs and point charges. The wave function of the ions are in parallel iterated until convergence in energy is achieved. This procedure is not variational, so only energy stabilisation is desired.
The hosts used in the optimisation of AIMPs with SCEPIC can be obtained from any program, but the embedding program described in the previous section is the recommended program to use.
Keywords related to the program itself are given in Table and for the electronic-structure calculations of the embedded ions in Table .
When running SCEPIC, the Molcas environment variables MOLCAS, and SCEPIC_SCRATCH must be set. (In previous versions of SCEPIC, WorkDir was used
insted of SCEPIC_SCRATCH).
SCEPIC can be used with either MOLCAS (8.4
) and OpenMolcas (v21.06).
It might work with earlier versions, but this is not guaranteed.
Using Julias threading system, SCEPIC can be executed in parallel, allowing the wave function of each ion in every iteration to be computed simultaneously. While it can be used with an MPI-parallelised version of Molcas, SCEPIC disables Molcas parallelism by default to avoid conflicts. The parallelism in SCEPIC is only over the number of ions, i.e., there is no point in allocating more threads than there are ions. To execute scepic in parallel, either "export JULIA_NUM_THREADS=number_of_ions" or run the Julia-script with "julia -t number_of_ions name_of_script.jl".
Remember to add the SCEPIC source directory to the Julia load path before running it as "export JULIA_LOAD_PATH=path_to_SCEPIC_source:$JULIA_LOAD_PATH".
Example of a simple script to generate a nanocluster for MgO. This example is also present in the test-directory. Note the order of "ion_dict" and "input" when calling the scepic function. Since all values in "input" have default values, it is also legal to run SCEPIC as "scepic(ion_dict)".
using scepic_module
Mg_s = [31837.96390000, 4638.58435000, 1036.10000000, 289.40599600, 93.21576980, 33.12376270, 12.43182320, 2.79558836, 0.93001709, 0.10672631, 0.04011496]
Mg_p = [99.72445620, 22.83590540, 6.89188917, 2.24703391, 0.71995582, 0.11289831]
ion1 = Dict("element" => "Mg",
"position" => [0.0, 0.0, 0.0],
"charge" => 2,
"s_basis" => Mg_s,
"p_basis" => Mg_p,
"hamiltonian" => "rx2c",
"nucleus" => "finite",
"SCF" => "RHF",
"ksdft" => "pbe",
"AIMPlabel" => "Mg.ECP.author.0s.0s.0e-update",
"AIMPembedding" => "./Mg.aimpfield",
"xfield" => "./Mg.xfield",
"external_density" => "./Mg.density")
O_s = [5779.29889000, 857.71286100, 193.81320000, 54.37489950, 17.35126050, 5.93118269, 1.04013399, 0.30991762]
O_p = [17.74190910, 3.86037124, 1.04716001, 0.27554862]
ion2 = Dict("element" => "O",
"position" => [0.0, 0.0, 0.0],
"charge" => -2,
"s_basis" => O_s,
"p_basis" => O_p,
"hamiltonian" => "rx2c",
"nucleus" => "finite",
"SCF" => "RHF",
"ksdft" => "pbe",
"AIMPlabel" => "O.ECP.author.0s.0s.0e-update",
"AIMPembedding" => "./O.aimpfield",
"xfield" => "./O.xfield")
ion_dict = [ion1, ion2]
input = Dict("maxiter" => 2)
scepic(ion_dict, input)
Example of how the file(s) specified by AIMPembedding might look. Labels specified here must be matched by the AIMPlabel-keywords.
Basis set
O.ECP.author.0s.0s.0e-update / AIMPLIB
pseudocharge
AAAAAA 0.000000000000000 2.139250000000000 2.139250000000000 Angstrom
...
AAAAEC -2.139250000000000 -2.139250000000000 -4.278500000000000 Angstrom
End of Basis
Basis set
Mg.ECP.author.0s.0s.0e-update / AIMPLIB
pseudocharge
AAAAED 2.139250000000000 0.000000000000000 0.000000000000000 Angstrom
...
AAAAIG -4.278500000000000 -2.139250000000000 -4.278500000000000 Angstrom
End of Basis
Example of an xfield-file that can be used by SCEPIC. Note that the coordinates has to be in Å
26840 angstrom
8.557000000000000 2.139250000000000 2.139250000000000 -2.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000
2.139250000000000 8.557000000000000 2.139250000000000 -2.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000
2.139250000000000 2.139250000000000 8.557000000000000 -2.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000
....
36.367249999999999 36.367249999999999 36.367249999999999 0.015625000000000 0.000000000000000 0.000000000000000 0.000000000000000
Perl script EPIC can be used to prefer the default settings for generation of AIMPs. The script requires one parameter - XYZ formatted file with some extensions. This XYZ file (based on the standard XYZ file format: first line contains the number of atoms, the second line is a comment line. The following lines contain element label (which can contain an underscore sign and a label, e.g. O_a). In addition to atomic coordinates, information about the translation vectors should follow:
#Origin 0 0 0 #Vector1 2.0 0 0 #Vector2 0 2.0 0 #Vector3 0 0 2.0For the atoms, for which AIMPs should be constructed, an additional information must be added: the formal atomic charge and orbitals (all occupied shells of the ion).
Mg 0 0 0 +2 sp
The script has a small configuration section, which can be edited.
To run the script, MOLCAS environment variable must be defined. As well as SCEPIC_ROOT and SCEPIC_SCRATCH.
Compatibility with older versions of SCEPIC: in August 2022, some minor changes has been made in the SCEPIC code, which do influence this interface script.
After the run of the script, several files are generated, including a bash script, named runme. If necessary, the user can edit intermediate files before running runme file.
If you already have computed AIMPs you can use embedding program to generate atomic charges and a template for your Molcas input. But it is also possible to use embedding to prepare an input for SCEPIC program (to generate the basis set file). The latter option is described at the next section.
First, we need to create a coordinate file, describing the unit cell and the lattice. The format of this file is an extension of ordinary XYZ files, where an additional information is added in the comment line. The comment line should contain keyword LATTICE, followed by 9 numbers. Note that the default unit in XYZ files is Angstroms.
The final result will contain three areas: the quantum part (QP), potential part (PP), and charges (CP). The QP is set by the coordinate file in XYZ format. The coordinates in the quantum part should match the positions of the lattice. PP and CP are defined in the driver file, which is described below.
An example of the driver file can be found in the test/embedding/in directory. It contains compulsory parts: using embedding_factory definition of unit cell (uc_file) and quantum part (qm_file) and the definition of all elements and formal charges in the following form: charge_dict = Dict("Mg" => 2, "O" => -2)
There are two alternative ways to setup the size of the CP: by radius (Charges_Radius) or by extension of the unit cells (Charges_Size). These keywords are exclusive. The Charges_Radius keyword requires one number. It can be complemented by keyword Charges_Center, which defines the center of the sphere.
If Charges_Size is defined instead, it requires one, three, or six numbers, defining the translation of the unit cell (uniform in all directions, symmetrical, or arbitrary).
The following keywords are optional. MaxMultipole keyword defines the highest Multiple moment, to be eliminated.
The setup of AIMP region is similar to set up of CP. Two alternative keywords can be used: AIMP_Radius or AIMP_Size.
The complete example of driver.jl script:
using embedding_factory
# define the coordinates of the system
#
# Unit cell file (it must contain in the comment line LATTICE followed by 9 numbers
uc_file = "./MgO_bulk.xyz"
# Quantum part (ordinary xyz file)
qm_file = "./MgO_clust.xyz"
#definition of elements and charges
charge_dict = Dict("Mg" => 2, "O" => -2)
# region of charges
Charges_Radius=20
Charges_Center=[0.0,0.0,0.0]
Charges_Size=[4, 4, 4]
# region of AIMPs
AIMP_Radius=5
AIMP_Size=[2, 2, 2]
#Multipoles to eliminate. Max=10, Default=4
MaxMultipole=3
#template for Spherical based cluster
input = Dict("elmoment" => MaxMultipole,
"unit_cell_file" => uc_file,
"qm_file" => qm_file,
"cluster_radius" => Charges_Radius,
"aimp_radius" => AIMP_Radius,
"reference" => Charges_Center,
"charge_dict" => charge_dict)
#template for Unit cell based cluster
input_UC = Dict("elmoment" => MaxMultipole,
"unit_cell_file" => uc_file,
"qm_file" => qm_file,
"cluster_size" => Charges_Size,
"aimp_radius" => AIMP_Size,
"reference" => Charges_Center,
"charge_dict" => charge_dict)
embedding(input)
#embedding(input_UC)
Some practical notes: the increase of the cost of calculations is insignificant, so try to use an ample layer of AIMPs and atomic charges.
For each atom in the system we need to prepare xfield file and ampfield file. In order to do that we have to run individual calculations for each atom. The lattice coordinate file is the same as in previous example, but quantum part contains only one single atom. In case of Mg, since it is located at (0,0,0), the only modification we need is to change the coordinate file, keeping only one atom in it. To avoid confusion, you can copy the coordinate file, edit it and change the name of the file in the driver script.
In the driver file we have to set Charges_Radius and AIMP_Radius to be large enough. AIMP_Radius must me at least 6Å, and for the Charges - the electrostatic potential should converge.
run julia driver.jl command. Rename project.xfield to Mg.xfield. Rename project.inp to Mg.aimpfield and remove from the top of the file the first 4 lines (the basis set for the real atom).
Next, we need to repeat the procedure for O atom. The coordinate file is:
1 comment O 2.128242 0 0
Since O is not in the center of coordinates, we have to change the center of sphere in the driver script:
qm_file = "./O_clust.xyz" Charges_Radius=40 Charges_Center=[2.1282420,0.0,0.0] AIMP_Radius=10
run julia driver-emb.jl and as before rename project.xfield to O.xfield and project.inp to O.aimpfield (remove 4 first lines).
The simplest way to generate the driver script for scepic is to use EPIC. In the beginning of the script, there is definition for each ion, which include the following parameters: element name, position (should match the corresponding selection in embedding input, formal atomic charge, exponents for valence shells (s- exponents are compulsory, others are optional), relativistic hamiltonian, nucleus, indication of Restricted or Unrestricted Hartree-Fock hamiltonian, DFT functional, and links to the files, generated by embedding program.
The exponents can be found in the basis library file. For UHF hamiltonian, the Multiplicity must be set.
Example of the driver file for MgO
using scepic_module
Mg_s = [169973.87,24699.066,7175.4851,2478.7920,917.73259,351.75140,137.51557,54.430286,21.722569,8.7194484,3.5147166,1.4212132,.57607650,.23395250,.09515600,.03875100,.01550040]
Mg_p = [557.13078,135.48388,47.212958,18.000955,7.1146607,2.8630240,1.1640185,.47630470,.19573560,.08067830,.03332660,.01333060]
ion1 = Dict("element" => "Mg",
"position" => [ 0.00000000, 0.00000000, 0.00000000],
"charge" => +2,
"s_basis" =>Mg_s,
"p_basis" =>Mg_p,
"hamiltonian" => "rx2c",
"nucleus" => "finite",
"SCF" => "RHF",
"ksdft" => "PBE",
"AIMPlabel" => "Mg.ECP.author.0s.0s.0e-update",
"AIMPembedding" => "./Mg.aimpfield",
"xfield" => "./Mg.xfield")
O_s = [105374.95,15679.240,3534.5447,987.36516,315.97875,111.65428,42.699451,17.395596,7.4383090,3.2228620,1.2538770,.49515500,.19166500,.06708300]
O_p = [200.00000,46.533367,14.621809,5.3130640,2.1025250,.85022300,.33759700,.12889200,.04511200]
ion2 = Dict("element" => "O",
"position" => [ 2.12824200, 0.00000000, 0.00000000],
"charge" => -2,
"s_basis" =>O_s,
"p_basis" =>O_p,
"hamiltonian" => "rx2c",
"nucleus" => "finite",
"SCF" => "RHF",
"ksdft" => "PBE",
"AIMPlabel" => "O.ECP.author.0s.0s.0e-update",
"AIMPembedding" => "./O.aimpfield",
"xfield" => "./O.xfield")
ion_dict =[ion1, ion2]
input = Dict("maxiter" => 20)
scepic(ion_dict)#, input)
executing this script will generate ECP file.
This document was generated using the LaTeX2HTML translator Version 2021.2 (Released July 1, 2021)
The command line arguments were:
latex2html -split 0 main231127.tex
The translation was initiated on 2023-11-27