*PEQM¶
Polarizable embedding model
This input section controls calculations using the polarizable embedding (PE) model. In DIRAC it is possible to include the effects from a structured environment on a core molecular system using the polarizable embedding (PE) model. The current implemen- tation is a layered QM/MM-type embedding model capable of using advanced potentials that include an electrostatic component as well as an induction (polarization) component. The effects of the environment are included through effective operators that contain an em- bedding potential, which is a representation of the environment, thereby directly affecting the molecular properties of the core system. The wave function of the core system is opti- mized while taking into account the explicit electrostatic interactions and induction effects from the environment in a fully self-consistent manner. The electrostatic and induction components are modeled using Cartesian multipole moments and anisotropic dipole-dipole polarizabilities, respectively. The electrostatic part models the permanent charge distribu- tion of the environment and will polarize the core system, while the induction part also allows polarization of the environment. The environment response is included in the effec- tive operator as induced dipoles, which arise due to the electric fields from the electrons and nuclei in the core system as well as from the environment itself. It is therefore necessary to recalculate the induced dipoles according to the changes in the electron density of the core system as they occur in a wave function optimization. Furthermore, since the induced dipoles are coupled through the electric fields, it is necessary to solve a set of coupled linear equations. This can be done using either an iterative or a direct solver. This also means that we include many-body effects of the total system.
The multipoles and polarizabilities can be obtained in many different ways. It is possible to use the molecular properties, however, usually distributed/localized properties are used because of the better convergence of the multipole expansion. These are typically centered on all atomic sites in the environment (and sometimes also bond-midpoints), how- ever, the implementation is general in this sense so they can be placed anywhere in space. Currently, the PE library supports multipole moments up to fifth order and anisotropic dipole-dipole polarizabilities are supported. For multipoles up to and including third order (octopoles) the trace will be removed if present. Note, that the fourth and fifth order multi- pole moments are expected to be traceless. In case polarizabilities are included it might be necessary to use an exclusion list to ensure that only relevant sites can polarize each other. The format of the POTENTIAL.INP file is demonstrated below.
Keywords¶
.POTENTIAL¶
This option can be used to specify a non-default name of the POTENTIAL.INP in- put file that contains the embedding potential parameters, e.g.:
.POTENT
peqm.pot
The default name is POTENTIAL.INP.
.DIRECT¶
Use the direct solver to determine the induced dipole moments. It will explicitly build a classical response matrix of size \(3S(3S + 1)/2\), where \(S\) is the number of polarizable sites and is therefore only recommendable for smaller molecular systems. Note that this solver is not parallelized. The default is to use the iterative solver
.ITERATIVE¶
Use the iterative solver to determine the induced dipole moments. This is the default. The convergence threshold defaults to \(1.0\cdot 10^{-8}>\sum^S_{s=1}|\mu^{(k)}_s-\mu^{(k-1)}_s|\), where \(k\) is the current iteration, but can also be provided with this option:
READ (LUCMD, *) THRITER (optional)
.BORDER¶
Controls the handling of the border between the core molecular system and its en- vironment described by the embedding potential.:
1) READ (LUCMD, *) BORDER_TYPE, RMIN, AUORAA
2) READ (LUCMD, *) BORDER_TYPE, REDIST_ORDER, RMIN, AUORAA, NREDIST
- There are two mutually exclusive schemes:
BORDER_TYPE = REMOVE
BORDER_TYPE = REDIST
The first option will remove all multipoles and polarizabilities that are within the given distance “RMIN” from any atom in the core molecular system. The “AUORAA” variable specifies whether the distance threshold is given in ångström (“AA”) or bohr (“AU”). The second option will redistribute parameters that are within the given threshold “RMIN” from any atom in the core system to nearest sites in the environment. The order of multipoles up to which will be redistributed is determined by the “REDIST_ORDER” variable, e.g. “REDIST_ORDER = 1” means that only charges will be redistributed and all other pa- rameters removed, “REDIST_ORDER = 2” means charges and dipoles are redistributed and so on. The sign of “REDIST_ORDER” specifies if the polarizabilities are redistributed. Positive means that the polarizabilities are removed and negative means redistributed. The number of sites that parameters on a given site are redistributed to is determined by the “NREDIST” variable which can be between 1 and 3. The default is to redistribute charges within 2.2 bohr to its nearest site and removing all other parameters.
.DAMP INDUCED¶
Damp the electric fields from induced dipole moments using Thole’s exponential damp- ing scheme.:
READ (LUCMD, *) IND_DAMP (optional)
The default damping coefficient is the standard 2.1304.
.DAMP MULTIPOLE¶
Damp the electric fields from permanent multipole moments using Thole’s expo- nential damping scheme.:
READ (LUCMD, *) MUL_DAMP (optional)
This option requires polarizabilities on all sites with permanent multipole moments. The default damping coefficient is the standard 2.1304.
.DAMP CORE¶
Damp the electric fields from the electrons and nuclei in the core region based on Thole’s exponential damping scheme:
READ (LUCMD, *) CORE_DAMP (optional)
READ (LUCMD, *) NALPHAS
DO I = 1, NALPHAS
READ (LUCMD, *) ISO_ALPHA(I)
END DO
The damping coefficient is optional unless isotropic polarizabilities are present. The default damping coefficient is the standard 2.1304. Standard polarizabilities from [vanduijnen1998] have been implemented and will be used if none are given as input. However, only H, C, N, O, F, S, Cl, Br and I are available.
.GSPOL¶
Activate the ground-state polarization approximation, i.e. freeze the embedding po- tential according to the ground-state density. This means that the polarizable envi- ronment is self-consistently relaxed during the optimization of the ground-state den- sity/wave function of the core molecular system and then kept frozen in any following response calculations
.NOMB¶
Remove many-body effects in the environment. This is done by deactivating interac- tions between inducible dipole moments.
.RESTART¶
Use any existing files to restart calculation.
.CUBE¶
Create cube file for the core molecular system containing the electrostatic potential due to the final converged polarizable embedding potential.:
1) READ (LUCMD, *) STD_GRID
2) READ (LUCMD, *) OPTION
2) IF (OPTION == ’GRID’) THEN
2) READ (LUCMD, *) XSIZE, XGRID, YSIZE, YGRID, ZSIZE, ZGRID
2) END IF
IF (OPTION == ’FIELD’) THEN
FIELD = .TRUE.
END IF
The grid density can be spec- ified either using 1) standard grids (COARSE (3 points/bohr), MEDIUM (6 points/bohr) or FINE (12 points/bohr) and in all cases 8.0 bohr are added in each direction) or 2) the GRID option which gives full control of the cube size and density, and requires an additional input line specifying the extent added (in bohr) in each direction and density (in points/bohr). The default grid density is MEDIUM and default cube size is the extent of the molecule plus 8.0 bohr in plus and minus each Cartesian coordi- nate. If the FIELD option is given then also three cube files containing the Cartesian components of the electric field from the embedding potential will be created.
.ISOPOL¶
Converts all anisotropic polarizabilities into isotropic polarizabilities.
.ZEROPOL¶
Remove all polarizabilities
.ZEROMUL¶
Remove all multipoles of order ZEROMUL_ORDER and up.:
READ (LUCMD, *) ZEROMUL_ORDER (optional)
Remove all multipoles of order ZEROMUL_ORDER and up. The default is to remove all dipoles and higher-order multipoles (i.e. ZEROMUL_ORDER = 1).
.VERBOSE¶
Verbose output. Currently this will print the final converged induced dipole moments.
.DEBUG¶
Debug output. Prints the total electric field and the induced dipole moments in each iteration. WARNING: for large systems this will produce very large output files.
The potential input format¶
The POTENTIAL.INP file is split into three sections: @COORDINATES, @MULTIPOLES and @POLARIZABILITIES. The format is perhaps best illustrated using an example:
! potential: 3_h2o.pot
! ------------------
! Multipole moments:
! ------------------
! TYPE: DALTON-LoProp
! METHOD: DFT
! XC-FUN: CAMB3LYP
! BASIS SET: loprop-6-31+G*
! -----------------
! Polarizabilities:
! -----------------
! TYPE: DALTON-LoProp
! METHOD: DFT
! XC-FUN: CAMB3LYP
! BASIS SET: loprop-6-31+G*
! -----------------
@COORDINATES
9
AA
O 2.36500000 -0.26700000 1.17000000
H 2.26100000 -1.10300000 1.62300000
H 3.13200000 -0.39000000 0.61100000
O 0.60500000 -1.46500000 -2.08000000
H 1.56100000 -1.47200000 -2.03300000
H 0.33300000 -0.92000000 -1.34100000
O -2.76000000 -0.37900000 1.04500000
H -3.11300000 -0.71800000 0.22300000
H -3.08300000 -0.98500000 1.71200000
@MULTIPOLES
ORDER 0
9
1 -0.71214000
2 0.35613000
3 0.35601000
4 -0.71172000
5 0.35592000
6 0.35580000
7 -0.71195000
8 0.35605000
9 0.35589000
ORDER 1
9
1 0.18080000 -0.26163000 -0.02886000
2 0.02828000 0.15363000 -0.09058000
3 -0.14340000 0.01302000 0.10892000
4 0.18641000 0.14653000 0.21409000
5 -0.18048000 0.00772000 -0.00000000
6 0.06166000 -0.10108000 -0.13640000
7 -0.18419000 -0.25746000 -0.04234000
8 0.06167000 0.05571000 0.16046000
9 0.05571000 0.10834000 -0.13341000
ORDER 2
9
1 -4.26483000 0.09135000 -0.61915000 -4.18829000 -0.42704000 -4.30439000
2 -0.44339000 0.02690000 -0.02255000 -0.13348000 -0.16450000 -0.36591000
3 -0.19077000 -0.05251000 -0.18965000 -0.43313000 0.03146000 -0.31920000
4 -3.75473000 -0.26370000 -0.29217000 -4.64173000 0.46544000 -4.36137000
5 -0.04700000 0.00328000 0.02886000 -0.45150000 0.00785000 -0.44444000
6 -0.41521000 -0.05871000 -0.07885000 -0.32164000 0.18442000 -0.20640000
7 -4.75797000 0.31610000 0.08263000 -4.49207000 -0.18934000 -3.50706000
8 -0.39388000 0.06405000 0.12905000 -0.39012000 0.12442000 -0.15864000
9 -0.40305000 0.09738000 -0.09242000 -0.27989000 -0.17467000 -0.26011000
@POLARIZABILITIES
ORDER 1 1
9
1 4.27640000 0.35715000 0.52792000 3.89073000 0.26437000 4.59489000
2 0.42029000 0.22994000 0.06545000 2.24667000 -0.89940000 0.97487000
3 1.98744000 -0.26213000 -0.97060000 0.39136000 0.31466000 1.26588000
4 3.86307000 -0.04205000 -0.13378000 4.67372000 -0.64393000 4.23650000
5 2.71461000 0.04981000 0.19272000 0.51208000 -0.15905000 0.42003000
6 0.43352000 -0.33472000 -0.47532000 1.31851000 0.93618000 1.89744000
7 4.64952000 -0.67410000 -0.13594000 4.15079000 0.04738000 3.96576000
8 0.88018000 0.20475000 0.75848000 0.75651000 0.80933000 2.00627000
9 0.82484000 0.41178000 -0.61423000 1.44166000 -1.04482000 1.38018000
EXCLISTS
9 3
1 2 3
2 1 3
3 1 2
4 5 6
5 4 6
6 4 5
7 8 9
8 7 9
9 7 8
@COORDINATES¶
The coordinates section follows the standard XYZ file format so that the environment can be easily visualized using standard programs. The first line in gives the total number of sites in the environment and the second line specifies whether the coordinates are given in ångström (AA) or bohr (AU). The rest of the coordinates section is a list of the sites in the environment where each line contains the element symbol and x-, y- and z-coordinates of a site. If a site is not located on an atom, e.g. if it is a bond-midpoint, then the element symbol should be specified as X. The listing also gives an implicit numbering of the sites, so that the first line is site number one, the second line is site number two and so on. This numbering is important and used in the following sections.
@MULTIPOLES¶
The multipoles section is subdivided into the orders of the multipoles, i.e. ORDER 0 for monopoles/charges, ORDER 1 for dipoles and so on. For each order there is a number specifying the number of multipoles of that specific order. Note, that this number does not have to be equal to the total number of sites. This is followed by a list of multipoles where each line gives the multipole of a site. The lines begin with a number that specifies which site the multipole is placed. Only the symmetry-independent Cartesian multipoles (given in a.u.) should be provided using an ordering such that the components are stepped from the right, e.g. xx xy xz yy yz zz or xxx xxy xxz xyy xyz xzz yyy yyz yzz zzz. Note, that the multipoles should in general be traceless, however, for multipoles up to and including octopoles (i.e. ORDER 3) the trace is removed if present. Furthermore, the current implementation is limited to fifth order multipoles.
@POLARIZABILITIES¶
The polarizabilities section is also subdivided into orders, i.e. ORDER 1 1 for dipole-dipole polarizabilities, which is the only type supported in the current release. The format is the same as for multipoles, i.e. first line contains number of polarizabilities which is followed by a list of the polarizabilities using the same ordering as the multipoles. The polarizabilities should also be given in a.u. In addition, there is also the exclusion lists (EXCLISTS section). Here the first line gives the number of lists (i.e. the number of lines) and the length of the exclusion lists (i.e. the number of entries per line). The exclusion lists specify the polar- ization rules. There is a list attached to each polarizable site that specifies which sites are not allowed to polarize it, e.g. 1 2 3 4 5 means that site number 1 cannot be polarized by sites 2, 3, 4 and 5.