Size: 50569
Comment:
|
Size: 50457
Comment:
|
Deletions are marked like this. | Additions are marked like this. |
Line 1: | Line 1: |
## page was renamed from AJMGrpOnly/camcasp/pols-cn ## page was renamed from ajm/camcasp/pols-cn |
|
Line 50: | Line 48: |
To calculate polarization energies we need the distributed multipole moments in addition to the distributed polarizabilities. Additionally, we need to determine the damping. For this see the sections on [[ajm/camcasp/polarization|Polarization]] and [[ajm/camcasp/charge-transfer|Charge-Transfer]]. | To calculate polarization energies we need the distributed multipole moments in addition to the distributed polarizabilities. Additionally, we need to determine the damping. For this see the sections on [[AJMPublic/camcasp/polarization|Polarization]] and [[AJMPublic/camcasp/charge-transfer|Charge-Transfer]]. |
Line 64: | Line 62: |
$$ E^{(2)}_{\rm disp}({\rm DM}) = - \sum_{a,b>a} \sum_{n=6}^{N_{\rm max}^{ab}} f_{n}(\beta^{ab} r_{ab}) \frac{C_{n}^{ab}(\Omega_{ab})}{r_{ab}^{n}} $$ |
$$E^{(2)}_{\rm disp}({\rm DM}) = - \sum_{a,b>a} \sum_{n=6}^{N_{\rm max}^{ab}} f_{n}(\beta^{ab} r_{ab}) \frac{C_{n}^{ab}(\Omega_{ab})}{r_{ab}^{n}}$$ |
Line 71: | Line 68: |
\begin{equation} f_{n}(\beta^{ab} r_{ab}) = 1 - \exp(-\beta^{ab} r_{ab}) \sum_{k=0}^{n}\frac{(\beta^{ab} r_{ab})^{k}}{k!} \end{equation} where $\beta^{ab}$ is the damping constant which may be angle-dependent, but is often assumed to be a simple constant dependent on the pair of sites $(a,b)$. |
$\begin{equation} f_{n}(\beta^{ab} r_{ab}) = 1 - \exp(-\beta^{ab} r_{ab}) \sum_{k=0}^{n}\frac{(\beta^{ab} r_{ab})^{k}}{k!} \end{equation}$ where $\beta^{ab}$ is the damping constant which may be angle-dependent, but is often assumed to be a simple constant dependent on the pair of sites $(a,b)$. |
Line 78: | Line 74: |
$$ \beta^{AB}_{\rm disp} = \sqrt{2 I_{\rm A}} + \sqrt{2 I_{\rm B}} $$ |
$$ \beta^{AB}_{\rm disp} = \sqrt{2 I_{\rm A}} + \sqrt{2 I_{\rm B}}$$ |
Contents
Navigation:
Distributed molecular properties
Introduction
CamCASP makes it very easy to calculate molecular properties like:
- Distributed multipole moments
- Distributed (local) polarizabilities, and
- Distributed (two-centred) dispersion coefficients
The multipoles are usually calculated using Anthony Stone's GDMA2 method. The GDMA code is strongly interfaced with CamCASP and is able to use wavefunctions read in via the CamCASP interfaces to various SCF/DFT codes like GAMESS(US), NWChem, DALTON2 and Gaussian (this is still a partial interface only).
Localized distributed polarizabilities are calculated using the Williams-Stone-Misquitta algorithm, and from these, the two-centerd dispersion models can be calculated.
As part of the polarizability calculation CamCASP calculates non-local (distributed) polarizabilities, and these can be use directly to calculate dispersion energies using the Dispersion programme. While this approach is very accurate and includes charge-flow effects that are important in low-dimensional delocalized systems (1-D wires), it is also quite computationally expensive to use these non-local models (the dispersion energy involves a quadruple sum over sites!). We will describe this sort of calculation elsewhere. For now we focus on the creation of local models.
Warning
CamCASP is in a state of change (goes without saying). This tutorial was written for version 5.9.x of the code, but from 6.0.x the format of the command files has changed - for the better.
Consequently, I have included both sets of files in the examples below. Please ensure that you are using files consistent with your version of the code!
Papers and Theory
I'll put in the relevant theoretical details later; for now, here are the papers that decsribe the methods:
Distributed Multipole Analysis: Stability for Large Basis Sets, Stone, J. Chem. Theory Comput., 1, 1128–1132 (2005). This is the GDMA2 method for distributed multipoles.
Distributed polarizabilities obtained using a constrained density-fitting algorithm, Misquitta and Stone, J. Chem. Phys. 124, 024111 (2006). These non-local polarizabilities form the basis for the WSM polarizabilities.
Accurate Induction Energies for Small Organic Molecules: 1. Theory, Misquitta and Stone, J. Chem. Theory Comput., 4, 7-18 (2008). Describes the Williams-Stone-Misquitta (WSM) polarizability models
Accurate Induction Energies for Small Organic Molecules. 2. Development and Testing of Distributed Polarizability Models against SAPT(DFT) Energies//, Misquitta, Stone and Price, J. Chem. Theory Comput., 4, 19-32 (2008). Examples of the WSM polarizabilities.
Dispersion energies for small organic molecules: first row atoms, Misquitta and Stone, Mol. Phys., 106, 1631-1643 (2008). Here we use the WSM models to create dispersion models.
Goals
Multipole models
These will be of the form: $Q^{a}_{lm}$, where $a$ is a site and $lm$ are the angular momentum labels. These are normally defined in the laboratory frame. They can be rotated into a local frame using the ORIENT program, but this is not normally done.
WSM polarizabilities
These are of the form : $\alpha^{a}_{lm,l'm'}$, where $a$ is a site and $lm$/$l'm'$ are the angular momentum labels. Since these are local polarizabilities, they depend on single sites only. These are normally defined in a local axis frame that is defined before the localization step in the WSM procedure. The axis frame should reflect the local symmetry, typically with the $z$-axis defined to be along terminal bonds.
To calculate polarization energies we need the distributed multipole moments in addition to the distributed polarizabilities. Additionally, we need to determine the damping. For this see the sections on Polarization and Charge-Transfer.
The WSM polarizabilities are always an approximation. It is the non-local polarizabilities that provide the full polarizability description. In localizing these we necessarily make an approximation. For many systems (those with large HOMO-LUMO gaps) the approximation may be a very good one, but for others (low dimensional systems with small gaps) the true physical polarizability description may indeed be non-local. In the latter case, the WSM description will fail.
However, when the WSM localization works, the WSM method gives us, in a well-defined sense, the best polarizability description given the limitations of the model. This is important. So, if we wish to have only rank 1 terms (dipole-dipole polarizabilities) on all atoms, then the WSM rank 1 model will be the best rank 1 model possible in the sense that it would be best able to reproduce the point-to-point polarizabilities. Likewise, if an isotropic polarizability description was needed, we'd best the most accurate model possible.
Important
The WSM models do not favour any particular type of interaction or orientation: The models will seek to best reproduce all the point-to-point polarizabilities which are uniformly distributed on a pseudo-random set of points around the molecule. The only way to bias the model is to bias the choice of point-to-point polarizabilities. So, for very specific application, when a model is desired to work for very specific interactions - say to favour just the hydrogen-bonded, or stacking orientations - the WSM models may not be the most appropriate.
WSM dispersion models
The dispersion models we need are in the form (using the SAPT notation; here the distributed multipole expansion is denoted by //DM//):
$$E^{(2)}_{\rm disp}({\rm DM}) = - \sum_{a,b>a} \sum_{n=6}^{N_{\rm max}^{ab}} f_{n}(\beta^{ab} r_{ab}) \frac{C_{n}^{ab}(\Omega_{ab})}{r_{ab}^{n}}$$
where $a$ and $b$ label sites, $f_{n}$ is a damping function of order $n$, $r_{ab}$ is the inter-site distance, $\Omega_{ab}$ are angular variables describing the orientation of the local axes on sites $a$ and $b$ relative the distance vector joining them, $N_{\rm max}^{ab}$ is the maximum order of the expansion for the particular site pair, and $C_{n}^{ab}(\Omega_{ab})$ are the angular-dependent dispersion coefficients.
We will typically use the Tang-Toennies incomplete Gamma functions for the damping:
$\begin{equation} f_{n}(\beta^{ab} r_{ab}) = 1 - \exp(-\beta^{ab} r_{ab}) \sum_{k=0}^{n}\frac{(\beta^{ab} r_{ab})^{k}}{k!} \end{equation}$
where $\beta^{ab}$ is the damping constant which may be angle-dependent, but is often assumed to be a simple constant dependent on the pair of sites $(a,b)$.
Often, even the site dependence is dropped and $\beta^{ab}$ is assumed to depend on the molecular types alone. Based on a simple idea, we have demonstrated that a good choice for the damping parameter for the //dispersion expansion only// is to use
$$ \beta^{AB}_{\rm disp} = \sqrt{2 I_{\rm A}} + \sqrt{2 I_{\rm B}}$$
where $I_{\rm A}$ and $I_{\rm B}$ are the vertical (first) ionization energies for the interacting molecules.
Our goal here is to use the WSM frequency-denpendent polarizabilities to calculate the $C_{n}^{ab}(\Omega_{ab})$ dispersion coefficients.
Important
We have recently started using damping models based on the ISA atomic ionization energies (those estimated from the ISA atomic domains). More on this later.
Distributed molecular properties with CamCASP 5.9.x and earlier
Warning
This section is for CamCASP 5.9.x and earlier versions only!
Files and Quantities
You will need:
- Molecular geometry in either XYZ format or the Z-matrix format.
- This should be a closed-shell system!!!
For calculations with an asymptotically corrected XC functional you will need the first vertical IP. Either use an experimental value or calculate it using the PBE0 functional in a $\Delta\mbox{DFT}$ procedure.
Files needed:
- Cluster file for the properties calculation.
- Axis file for local axis definitions (if not used, all properties will be calculated in the global (lab) axis system).
We will illustrate this with the water molecule. Here is the cluster file for water (file name: H2O.clt):
Title H2O properties calculation Global Units Bohr Degrees Overwrite Yes End Molecule H2O I.P. 0.4638 a.u. ( NIST Chem Webbook ) ! Vibrationally averaged geometry from ! Mas and Szalewicz [J Chem Phys 104, 7606 (1996)]. Units Bohr O 8.0 0.0000000000 0.0000000000 0.0000000000 Type O H1 1.0 -1.4536519600 0.0000000000 -1.1216873200 Type H H2 1.0 1.4536519600 0.0000000000 -1.1216873200 Type H End Run-Type Molecule H2O File-prefix H2O_aTZ ! Basis aVTZ Aux-Basis aVTZ ! Functional PBE0 AC MULTPOLE Kernel ALDA+CHF SCF-Code DALTON2006 ( As I still use the old DALTON ) ! Orient files for localization and display Process file Sites file Interface file End Finish
Some of the options are defaults (PBE0, DALTON) and some would be automatically set given others (the aug-cc-pVTZ auxiliary basis would automatically be set with the aug-cc-pVTZ main basis). See the CamCASP User's Guide for details of the Cluster input.
The axis file takes the form (file name H2O.axes):
Axes O z between H1 and H2 y from H1 to H2 H1 z from H1 to O y from H1 to H2 H2 z from H2 to O y from H2 to H1 End
See the ORIENT manual for detail of the syntax.
This job finishes rather quickly. The results will be in a directory //H2O_aTZ/OUT/ //. Have a look at the CamCASP output in file //H2O_aTZ/OUT/H2O_aTZ_camcasp.out//. The total polarizability tensor should be something like:
Order: 1 by 1 9.391882 0.000000 0.000000 0.000000 10.034214 0.000000 0.000000 0.000000 8.750848 Isotropic polarizability: 9.39231483 Anisotropic polarizability: 1.11142748
The formatting will not be the same as above.
The //OUT/ // directory will contain a few more files:
$ ls H2O_aTZ/OUT/ H2O_aTZ-data-summary.data H2O_aTZ.out H2O_aTZ_0.0005_1000_f11_NL4.pol H2O_aTZ_DMA2_L4.mom H2O_aTZ_camcasp.out H2O_aTZ_dalton.out H2O_daTZ_lim2.0_4.0_p2000_f11.p2p
The non-local polarizabilities are in file //H2O_aTZ_0.0005_1000_f11_NL4.pol//, the distributed multipole moments in //H2O_aTZ_DMA2_L4.mom// and the point-to-point polarizabilities in file //H2O_daTZ_lim2.0_4.0_p2000_f11.p2p//.
Localizing the non-local polarizabilities
In this step we will localize the non-local polarizabilities and create a dispersion model for the water dimer. We will create a rank 3 model. By default the scripts will limit polarizabilities on hydrogen atoms to rank 1. If you wish to increase this use the //--hlimit// option to the localize script.
$ mkdir L3 $ cd L3 $ cp ../H2O.clt . $ cluster < H2O.clt ( we need the various interface files Cluster creates) $ ln -s ../H2O_aTZ/OUT . ( we need to make this link as the localize script need to find the point-to-point and the non-local polarizabilities ) $ localize.py --limit 3 --axes ../H2O.axes H2O_aTZ Using polarizability file OUT/H2O_aTZ_0.0005_1000_f11_NL4.pol Using p2p file OUT/H2O_aTZ_lim2.0_4.0_p2000_f11.p2p Splitting OUT/H2O_aTZ_0.0005_1000_f11_NL4.pol into individual frequencies .../usr/bin/csplit --prefix=H2O_aTZ_NL4_ --suffix-format=%03d.pol --elide-empty-files --quiet OUT/H2O_aTZ_0.0005_1000_f11_NL4.pol /# INDEX/ '{*}' done Localizing polarizabilities using LW procedure ... 000 001 002 003 004 005 006 007 008 009 010 ... done Splitting the point-to-point polarizability file ... Splitting point-to-point pols for H2O_aTZ Frequencies : 0 to 10 Number of points in grid : 2000 Total p2p pols per frequency : 2001000 Split p2p pols in files with prefix : H2O_aTZ_010 done Refining the local polarizabilities 000 001 002 003 004 005 006 007 008 009 010 ... done Calculating the dispersion coefficients ... ...done. Dispersion coefficients are in H2O_aTZ_ref_wt4_L3_casimir.out The dispersion potential, in Orient form, is in H2O_aTZ_wt4_L3_C12.pot $
The localize.py script does quite a number of things. In brief:
- Split the point-to-point polarizability file into files for each of the 11 frequencies (Static + 10 imaginary frequencies. The latter for the integration to get the dispersion coefficients using the Casimir-Polder integrals).
- Run Cluster to generate the ORIENT input files for localization.
- Run ORIENT to localize the non-local polarizabilities (using, in this example and by default, the Lillestolen-Wheatley localization scheme).
- Run the Process code to write out the PFIT commands for the refinement stage.
- Run PFIT to refine the localized polarizabilities.
- Run Process to write out the Casimir commands.
- Run Casimir to calculate the dispersion coefficients.
And the quantities we want are:
- The localized and refined WSM rank 3 polarizabilities in file //H2O_daTZ_ref_wt4_L3_static.pol//
- and the dispersion model in file //H2O_daTZ_wt4_L3_C12.pot//
Important
The localize.py script attempts to do a lot of steps and things can go wrong. If something does, use the //--keep// option to get it to keep all intermediate files so you get an idea of when and where it failed. Start by looking into the results of the localization and refinement steps. The latter will be in files of the sort //H2O_daTZ_ref_wt4_L3_000.out//. The total polarizabilities are printed out at the end. Check to see if they make sense.
Warning
Make sure your local axis definition reflects the symmetry of the system. If it doesn't, the total molecular polarizability will come out wrong. As yet there isn't a fool-proof way of doing this. Best to draw a diagram and make sure your axis system satisfies the point group (full or sub-group) of the system.
For a rank 1 isotropic model we would use (in this case we use a d-aug-cc-pVTZ basis):
$ localize.py --axes H2O.axes --limit 1 H2O_daTZ
This gives the following WSM static polarizability model (done with a d-aug-cc-pVTZ basis):
$ more H2O_daTZ_ref_wt4_L1_static.pol # Static polarizabilities O O 0.000 0.0000000 0.0000000 0.0000000 0.000 6.6033732 0.0000000 0.0000000 0.000 0.0000000 6.2713511 0.0000000 0.000 0.0000000 0.0000000 6.5586589 H1 H1 0.000 0.0000000 0.0000000 0.0000000 0.000 2.2630552 0.42434356E-01 0.0000000 0.000 0.42434356E-01 1.0983324 0.0000000 0.000 0.0000000 0.0000000 1.3045396 H2 H2 0.000 0.0000000 0.0000000 0.0000000 0.000 2.2630552 0.42434356E-01 0.0000000 0.000 0.42434356E-01 1.0983324 0.0000000 0.000 0.0000000 0.0000000 1.3045396
For an //isotropic// rank 3 model use
$ localize.py --limit 3 --isotropic H2O_daTZ
In the latter case, the dispersion model is
$ more H2O_daTZ_wt4_L3iso_C12iso.pot O O C6 C7 C8 C9 C10 C11 C12 00 00 0 24.34135 0.0 537.4367 0.0 15579.54 0.0 332073.2 End H O C6 C7 C8 C9 C10 C11 C12 00 00 0 4.381638 0.0 48.00146 0.0 900.4451 End H H C6 C7 C8 C9 C10 C11 C12 00 00 0 0.8002536 End
The //00 00 0// entries are angular momenta labels, all zero for static coefficients. Notice that the odd terms are zero for an isotropic dispersion model.
Important
All results are in atomic units.
Important
While isotropic dispersion models work reasonably well (see Misquitta and Stone, 2008), it may be better to use an anisotropic static polarizability model. Also, if you wish to use a simple $C_6$ dispersion model, consider scaling it to mimic the effects of the missing higher-ranking terms. This is an approximation, but it has been shown to give reasonable result in our 2008 Mol. Phys. paper.
Distributed molecular properties with CamCASP 6.0.x and later
Warning
This section is for CamCASP 6.0.x
Files and Quantities
We perform the properties calculation in two stages:
- First perform an ISA calculation, then
- perform an ISA-Pol calculation.
This can be done in one step, but as we are still learning about the ISA, I find it best to not treat the ISA as a black-box. Instead I always look at the ISA solution to see if it makes sense.
You will need:
- Molecular geometry in either XYZ format or the Z-matrix format.
- This should be a closed-shell system!!!
For calculations with an asymptotically corrected XC functional you will need the first vertical ionization potential (IP), and possibly the AC shift, $\Delta = {\rm IP} + \epsilon_{\rm HOMO}$. For the IP, either use an experimental value or calculate it using the PBE0 functional in a $\Delta\mbox{DFT}$ procedure. The energy of the HOMO level should be calculated with your functional of choice and without the AC!
Files needed:
- Cluster file for the properties calculation.
- Axis file for local axis definitions (if not used, all properties will be calculated in the global (lab) axis system).
We will illustrate this with the water molecule. Here is the cluster file for water (file name: H2O.clt):
Title H2O : ISA calculation with CamCASP 6.0.x Global Units Bohr Degrees Overwrite Yes End Molecule H2O I.P. 0.4638 a.u. ( NIST Chem Webbook ) AC-Shift 0.1311 a.u. ( for PBE0 with an aTZ basis ) ! Vibrationally averaged geometry from ! Mas and Szalewicz [J Chem Phys 104, 7606 (1996)]. Units Bohr O 8.0 0.0000000000 0.0000000000 0.0000000000 Type O H1 1.0 -1.4536519600 0.0000000000 -1.1216873200 Type H H2 1.0 1.4536519600 0.0000000000 -1.1216873200 Type H End Show H2O in XYZ format ( you could also use PDB as the format ) Run-Type Properties Molecule H2O Main-Basis aVTZ Type MC Aux-Basis aVTZ Type MC Cartesian Use-ISA-Basis ( this basis *could* be Spherical ) AtomAux-Basis aVQZ Type MC Spherical Use-ISA-Basis ( this *must* be Spherical ) ISA-Basis set2 Min-S-exp-H = 0.2 ( convergence tends to be better with this setting ) Func PBE0 ( other functionals possible, but this is probably the best ) AC CS00 ( [TH | LINEAR | GRAC] [CS00 | MULTIPOLE | LB94] : CS00 with NWChem) Kernel ALDA+CHF ( not needed for an ISA calculation ) SCF-code NWChem ( DALTON2015, NWCHEM, etc ) File-Prefix H2O_aTZ ( change this to something reasonable ) #METHOD isa-A ( this will use the file isa-A from the methods/ directory ) End Finish
The axis file takes the form (file name H2O.axes):
Axes O z between H1 and H2 y from H1 to H2 H1 z from H1 to O y from H1 to H2 H2 z from H2 to O y from H2 to H1 End
See the ORIENT manual for detail of the syntax.
Running CamCASP
Step (1): ISA
In this step we perform the DFT calculation, obtain the molecular orbitals, and calculate the ISA solution and ISA multipoles.
If the CamCASP code has been set up correctly, all you need to do is execute:
$ runcamcasp.py --clt H2O-ISA.clt -q bg --nproc 1 H2O_aTZ >& log &
where //H2O_aTZ// is the file prefix defined in the Cluster input file, //H2O-ISA.clt// is the name of the cluster input file (by default the scripts will look for a file with the name PREFIX.clt) and I have sent the job to the background using //-q bg//. For more details please see the User's Guide.
If all runs correctly, you will have a directory H2O_aTZ/OUT/ in which you should see the files There may be other files too. The CamCASP output is in Important The ISA solution is in file H2O_atoms.ISA. This file will be needed for any subsequent restart of the ISA. We will need it for the next step.
Now we need to calculate the ISA-Pol non-local polarizabilities. This will be done as follows: Run CamCASP only using the 'runcamcasp' command (note, there is no In the future we will streamline these steps, but for now you do need to do some things by hand. Here are some more details: - Go to directory Warning The above commands assume you have a directory structure like this: If your directory structure differs from this, then correct the paths used in the link commands (ln -s) accordingly. For example, a common directory structure is In this case the link commands would be:
The Cluster file for this calculation is very similar to the one use for the ISA calculation. There are some changes to the The RUN-TYPE will be Change the METHOD to The complete Cluster file should look like this: The axis file we will use is Warning It is very important to make sure that the axis file correctly reflects the symmetry of the system. In the above example we have indicated that the sites H1 and H2 both have the same TYPE H. This means that we have to use a local axis system in which these types are indeed the same and the above one does this as this axis system obeys the $C_{2v}$ symmetry of H2O.
This time we do not need to run NWChem as we already have the molecular orbitals, so rather than use the Warning One version of the Cluster program needs the job name specified using --job as: $ cluster --job H2O_aTZ < H2O-ISApol.clt This will not take long. CamCASP will first restart the ISA from the solution in The important files are: Warning I used a slightly different calculation to generate these examples. So file prefixes are sometimes H2O and other times H2O_aTZ. Keep this in mind! You should have the following total polarizability tensor: These were calculated with a larger basis, but the aug-cc-pVTZ basis should result in something like this.
The non-local polarizabilities can be localized using the Let's have a look at both isotropic and anisotropic models.
Here we will create a rank 3 isotropic local model for water. Make sure you are in the This will work, but I prefer to place the commands in a file (easy to keep a record and make modifications later), and also to perform the localization in a sub-directory as that way things are kept in order. So instead of the above use the following: Important Place the above commands in a file. Here we have called it loc-iso-command.bash, but you could use any name. Make this file executable using the chmod command: Now you can execute the commands in this file using: Here are the options used in the above: --subdir L3iso-wt3-coeff1e-3 : do the calculations in sub directory When executed you will get the following output: Warning The localization script will print out the localization settings for your reference. These are also printed out in the dispersion model file (below). Notice that there seems to be an axis file present. There will be one listed even if you didn't provide any as localize.py will create an empty file if it finds none present. And here is the content of the potential file: One way of assessing if this is a good model is to have a look at the output files of the form Look for the following:
The only difference with anisotropic models is that you need to have the axis file present in the work directory ( Warning Make sure you have copied or linked the H2O.axes file to the ISA-Pol/ directory! Here is the localization command: Execute this as before... And the output will be in directory The potential file now looks like this: Once again we should assess the localization by looking at the file Notice that the residuals are much smaller this time. Also, notice that the mean and anisotropic dipole-dipole polarizabilities are very close to those from CamCASP (above). Warning GOTCHAs If your axis system is inconsistent with the TYPE symmetry declared in the Cluster file, then the residuals will be quite a lot larger and the solution will be garbage.
The localize.py script works for homogeneous dimers. To get dispersion coefficients for a heterogeneous dimer, say water and Cl$^-$, you will need to first follow the above procedure to calculate the WSM models for the homogeneous case for both molecules. We then need to create the Casimir command file for the dimer. For this we will use the Process code. Then run Casimir and get our dispersion model. In a sense, all we need to do is give Casimir the WSM frequency-dependent polarizabilities in the right format. Let's say we have calculated the WSM polarizabilities and now have isotropic rank 3 static polarizabilities and dispersion coefficients for H$_2$O and Cl$^-$. We will make use of a few files from these calculations: WSM polarizabilities at all frequencies for each of the monomers: These are in the files //H2O_aTZ_ref_wt4_L3_0f10.pol// and //Clminus_ref_wt4_L3_0f10.pol// (I used the file prefix //Clminus// for the Cl$^-$ calculation). What we need to do is use the data in these files as the input to the Casimir code which uses the localized polarizabilities to calculate the dispersion coefficients. Unfortunately, at present, this requires creating an input file to the Process program by hand. Create a file //Clminus_H2O_casimir.prss// that contains the following: All the information in this file is contained in the files //Clminus_casimir.prss// and //H2O_aTZ_casimir.prss//. All we have done here is to merge them and add one additional line. Now create the Casimir command file using: This file looks like: Important The path to CamCASP is used in the last line where it points to the location of the real Clebsch-Gordon coefficients which are in file //realcg//. If this path is not correct or if you do not have read privileges to the data in file //readcg// the next step will fail. Now run the Casimir program: The output contains the dispersion model. It will be at the end of this file. For this example we get: Important If anisotropic polarizabilities are used, the resulting dispersion model will be in the axis system used for those polarizabilities!
A quick word about these: molecular properties like polarizabilities and dispersion coefficients require quite large and diffuse basis sets before they are sufficiently converged. The //aug-cc-pVTZ// basis set used in the above calculations is not good enough, particularly for the anion Cl$^-$. It is better to use a much more diffuse basis. Only bear in mind that there may not be an optimized auxiliary basis available for some of the very diffuse basis sets like the //d-aug-cc-pVTZ//. Some degree of experimentation with basis sets may be required.
We will normally use linear-response time-dependent (LR-TDDFT) DFT with PBE0 and the hybrid ALDA+CHF response kernel. But for some systems, for example anions like Cl$^-$, PBE0 may not be adequate because of the self-interaction error. The asymptotic correction can correct for this in some cases, but when it fails, you may need to use LC-PBE0 and use IP-tuning to determine the range-separation parameter. However this leads to a problem: the response kernel in CamCASP cannot yet be used consistently for range-separated functionals. So in this case please contact us to see what your options are.
The dispersion models should be tested against calculated SAPT(DFT) dispersion energies (use the total dispersion energy $E^{(2)}_{\rm disp,tot} = E^{(2)}_{\rm disp,pol} + E^{(2)}_{\rm disp,exch}$ in the comparisons). To do this we need the damping parameter. A good choice for the dispersion models is the expression given above with $\beta$ a function of the vertical ionization energies of the interacting molecules. Alternatively, this coefficient can be obtained by fitting to the SAPT(DFT) total dispersion energies. Model energies can be evaluated using the ORIENT program for dimer geometries specified in the angle-axis format. See the Orient: Energy Calculations section for more details. AIM-isosurface_H1.grid H2O_atoms.ISA H2O_ISA.mom w-H2--fn1-CONV.dat
AIM-isosurface_H2.grid H2O.cks warnings.log w-O--exp-analysis.dat
AIM-isosurface_O.grid H2O.out w-H1--exp-analysis.dat w-O--fn1-CONV.dat
data-summary.data H2O-ISA.clt w-H1--fn1-CONV.dat
H2O_1.00E-03_iso.grid H2O_ISA-GRID.mom w-H2--exp-analysis.dat
Step (2) : ISA-Pol
.py)
H2O_aTZ/. This is the work directory for this job. - Make directory ISA-Pol and go there. - Copy the H2O.clt: cp ../H2O.clt ./H2O-ISApol.clt. We will edit it as shown below. - Make links: ln -s ../*.movecs . and ln -s ../OUT/*.ISA . - Check: Do the links actually point to the files? Use ls -al to check. - Edit the H2O-ISApol.clt file as shown next.
/ISA-Pol/ Cluster file for ISA-Pol
PROPERTIES Title H2O : ISA-Pol calculation with CamCASP 6.0.x
Global
Units Bohr Degrees
Overwrite Yes
End
Molecule H2O
I.P. 0.4638 a.u. ( NIST Chem Webbook )
AC-Shift 0.1311 a.u. ( for PBE0 with an aTZ basis )
! Vibrationally averaged geometry from
! Mas and Szalewicz [J Chem Phys 104, 7606 (1996)].
Units Bohr
O 8.0 0.0000000000 0.0000000000 0.0000000000 Type O
H1 1.0 -1.4536519600 0.0000000000 -1.1216873200 Type H
H2 1.0 1.4536519600 0.0000000000 -1.1216873200 Type H
End
Show H2O in XYZ format ( you could also use PDB as the format )
Run-Type
Properties
Molecule H2O
Main-Basis aVTZ Type MC
Aux-Basis aVTZ Type MC Cartesian No-ISA-Basis ( CHANGE MADE HERE )
AtomAux-Basis aVQZ Type MC Spherical Use-ISA-Basis
ISA-Basis set2 Min-S-exp-H = 0.2
Func PBE0 ( other functionals possible, but this is probably the best )
AC CS00 ( [TH | LINEAR | GRAC] [CS00 | MULTIPOLE | LB94] : CS00 with NWChem)
Kernel ALDA+CHF ( must be this for the PBE0 functional )
SCF-code NWChem ( DALTON2015, NWCHEM, etc )
File-Prefix H2O_aTZ ( change this to something reasonable )
#METHOD isa-pol-from-isa-A-restart ( CHANGE MADE HERE )
End
Finish
Axes
O z between H1 and H2 y from H1 to H2
H1 z from H1 to O y from H1 to H2
H2 z from H2 to O y from H2 to H1
End
Running the ISA-Pol calculation
$ cluster < H2O-ISApol.clt
$ runcamcasp -q bg H2O_aTZ
H2O_atoms_001.ISA H2O_ISA-GRID_f11_NL4_fmtA.pol warnings.log w-H2--fn1-CONV.dat
H2O_atoms.ISA H2O_ISA-GRID_f11_NL4_fmtB.pol w-H1--exp-analysis.dat w-O--exp-analysis.dat
H2O_aTZ-A-asc.movecs H2O_ISA-GRID.mom w-H1--fn1-CONV.dat w-O--fn1-CONV.dat
H2O_aTZ.out H2O_lim2.0_4.0_p2000_f11.p2p w-H2--exp-analysis.dat
Spherical tensor form
Origin is ( 0.0000, 0.0000, 0.0000) BOHR
Order: 0 by 0
0.0000000E+00
Order: 0 by 1
0.3077651E-07 0.0000000E+00 0.0000000E+00
Order: 0 by 2
0.6482435E-07 0.0000000E+00 0.0000000E+00 0.3091289E-06 0.0000000E+00
Order: 1 by 1
0.9542880E+01 0.6063761E-05 0.0000000E+00
0.6063696E-05 0.9985474E+01 0.0000000E+00
0.0000000E+00 0.0000000E+00 0.9129133E+01
Isotropic polarizability: 0.9552495E+01
Anisotropic polarizability: 0.7417533E+00
Localization of the ISA-Pol polarizabilities
Isotropic local models
localize.py --limit 3 --wsmlimit 3 --hlimit 3 H2O_aTZ
localize.py --limit 3 --wsmlimit 3 --hlimit 3 \
--isotropic --subdir L3iso-wt3-coeff1e-3 \
--loc LW --weight 3 --weightcoeff 1e-3 \
H2O_aTZ
L3iso-wt3-coeff1e-3/ $ ./loc-iso-command.bash
Localisation settings for H2O_aTZ
Sub-directory L3iso-wt3-coeff1e-3-new
Axes file: H2O.axes
Pol file format: NEW
Limit: 3
WSM-Limit: 3
H-Limit: 3
Isotropic?: isotropic
Model file: H2O.pdef
Pol Cutoff: 0.0001
Loc algorithm: LW
Weight: 3
Weight coeff: 0.001
SVD threshold: 0.0
NoRefine?: False
Using polarizability file OUT/H2O_ISA-GRID_f11_NL4_fmtB.pol
Splitting OUT/H2O_ISA-GRID_f11_NL4_fmtB.pol into individual frequencies ... done
Localizing polarizabilities using LW procedure ...
000 001 002 003 004 005 006 007 008 009 010 ... done
Preparing to refine the local polarizabilities
Using axis definition file H2O.axes
Refining the local polarizabilities
Creating new polarizability model definition file H2O.pdef
000 001 002 003 004 005 006 007 008 009 010
Refinement finished
Refined localized polarizabilities, static and at imaginary frequency,
are in H2O_ref_wt3_L3iso_0f10.pol
Calculating the dispersion coefficients ... done
Dispersion coefficients are in H2O_ref_wt3_L3iso_casimir.out.
The dispersion potential, in Orient form, is in H2O_ref_wt3_L3iso_C12iso.pot.
! Localisation settings for H2O
! Axes file: H2O.axes
! Pol file format: NEW
! Limit: 3
! WSM-Limit: 3
! H-Limit: 3
! Isotropic?: True
! Model file: H2O.pdef
! Pol Cutoff: 0.0001
! Loc algorithm: LW
! Weight: 3
! Weight coeff: 0.001
! SVD threshold: 0.0
! NoRefine?: False
!
! Dispersion coefficients for H2O (hartree bohr^n)
O O C6 C7 C8 C9 C10 C11 C12
00 00 0 24.56813 0.0 487.0036 0.0 12225.23 0.0 229547.1
End
H O C6 C7 C8 C9 C10 C11 C12
00 00 0 4.297156 0.0 65.58407 0.0 1317.060 0.0 17518.63
End
H H C6 C7 C8 C9 C10 C11 C12
00 00 0 0.7630528 0.0 8.072538 0.0 112.2174 0.0 1262.786
End
...
...
Parameter values:
1 6.78769194 O_1_iso_A
2 24.82298591 O_2_iso_A
3 217.75336965 O_3_iso_A
4 1.35712844 H1_1_iso_A
5 2.20851916 H1_2_iso_A
6 6.64427990 H1_3_iso_A
Residuals per cent of range
Maximum 0.00130009 6.951
Minimum -0.00122429 -6.545
R.m.s. 0.00006801 0.364
Sum of squared residuals = 0.00926
Parameter Fitted Anchor Difference Penalty
O_1_iso_A 6.78769 7.08513 -0.29744 1.72795E-06
O_2_iso_A 24.82299 29.95894 -5.13595 2.93566E-05
O_3_iso_A 217.75337 217.73979 0.01358 3.88949E-12
H1_1_iso_A 1.35713 1.23371 0.12342 6.03978E-06
H1_2_iso_A 2.20852 2.33266 -0.12415 2.39268E-06
H1_3_iso_A 6.64428 6.63712 0.00716 1.13848E-09
...
...
Total molecular polarizability
00 z x y 20 21c 21s 22c 22s
0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
0.00000 9.50195 0.00000 0.00000 -6.08910 0.00000 0.00000 0.00000 0.00000
0.00000 0.00000 9.50195 0.00000 0.00000 -5.27331 0.00000 0.00000 0.00000
0.00000 0.00000 0.00000 9.50195 0.00000 0.00000 -5.27331 0.00000 0.00000
0.00000 -6.08910 0.00000 0.00000 48.63565 0.00000 0.00000 -9.93419 0.00000
0.00000 0.00000 -5.27331 0.00000 0.00000 56.69164 0.00000 0.00000 0.00000
0.00000 0.00000 0.00000 -5.27331 0.00000 0.00000 39.48512 0.00000 0.00000
0.00000 0.00000 0.00000 0.00000 -9.93419 0.00000 0.00000 46.44655 0.00000
0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 46.44655
Mean polarizability = 9.50195
Polarizability anisotropy = 0.00000
Finished
Anisotropic local models
localize.py --limit 3 --wsmlimit 3 --hlimit 3 \
--axes H2O.axes --subdir L3-wt3-coeff1e-3 \
--loc LW --weight 3 --weightcoeff 1e-3 \
H2O
$ ./loc-aniso-command.bash
! Localisation settings for H2O
! Axes file: H2O.axes
! Pol file format: NEW
! Limit: 3
! WSM-Limit: 3
! H-Limit: 3
! Isotropic?: False
! Model file: H2O.pdef
! Pol Cutoff: 0.0001
! Loc algorithm: LW
! Weight: 3
! Weight coeff: 0.001
! SVD threshold: 0.0
! NoRefine?: False
!
! Dispersion coefficients for H2O (hartree bohr^n)
O O C6 C7 C8 C9 C10 C11 C12
00 00 0 25.61645 0.0 570.7627 0.0 13821.21 0.0 263833.7
00 10 1 0.0 12.23063 0.0 214.5696 0.0 4057.924
00 11c 1 0.0 0.5356526E-02 0.0 0.8741571E-01 0.0 1.567216
00 20 2 -0.1775111 0.0 2.132830 0.0 -165.2042 0.0 -5261.844
00 21c 2 0.0 0.0 0.4491132E-02 0.0 0.2977902 0.0 8.375013
00 22c 2 -0.3082694 0.0 -117.7108 0.0 -2372.298 0.0 -50530.35
00 30 3 0.0 -3.660322 0.0 -70.86989 0.0 -1286.275
00 31c 3 0.0 -0.9719078E-03 0.0 -0.2925074E-01 0.0 -0.5761495
00 32c 3 0.0 6.174164 0.0 139.9797 0.0 2632.120
00 33c 3 0.0 -0.3764183E-02 0.0 -0.7991419E-01 0.0 -1.477058
...
...
...
...
Residuals per cent of range
Maximum 0.00029275 1.565
Minimum -0.00017147 -0.917
R.m.s. 0.00000528 0.028
Sum of squared residuals = 0.00006
Parameter Fitted Anchor Difference Penalty
O_10_10_A 6.93856 6.96757 -0.02901 1.69829E-08
O_10_20_A 0.28301 0.28126 0.00175 2.85077E-09
O_10_22c_A 1.29707 1.30912 -0.01205 5.35476E-08
O_10_30_A -1.92469 -1.92309 -0.00160 5.47390E-10
O_10_31c_A 0.00117 0.00058 0.00059 3.45231E-10
O_10_32c_A -5.06593 -5.12840 0.06247 1.42943E-07
O_10_33c_A 0.00030 0.00025 0.00005 2.40869E-12
O_11c_11c_A 6.75861 6.77750 -0.01889 7.60090E-09
O_11c_21c_A 2.60901 2.60524 0.00378 1.83090E-09
O_11c_31c_A 5.86913 5.86684 0.00229 1.47921E-10
...
...
there are lots of parameters this time...
...
...
Total molecular polarizability
00 z x y 20 21c 21s 22c 22s
0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
0.00000 9.52357 0.00000 0.00000 -4.87982 0.00000 0.00000 -2.84111 0.00000
0.00000 0.00000 9.97845 0.00000 0.00000 -10.02021 0.00000 0.00000 0.00000
0.00000 0.00000 0.00000 9.09659 0.00000 0.00000 -4.55987 0.00000 0.00236
0.00000 -4.87982 0.00000 0.00000 45.20779 -0.00055 0.00000 2.24896 0.00000
0.00000 0.00000 -10.02021 0.00000 -0.00055 57.65207 0.00000 0.00034 0.00000
0.00000 0.00000 0.00000 -4.55987 0.00000 0.00000 46.54062 0.00000 -0.00051
0.00000 -2.84111 0.00000 0.00000 2.24896 0.00034 0.00000 47.39755 0.00000
0.00000 0.00000 0.00000 0.00236 0.00000 0.00000 -0.00051 0.00000 44.94860
Mean polarizability = 9.53287
Polarizability anisotropy = 0.76384
Finished
Hetero-dimers : Dispersion models
TITLE PROCESS file to write the CASIMIR input for Cl- and H2O
Set Global-data
! This should be the PATH to your CamCASP installation:
CamCASP /home/alston/CamCASP/current
Units BOHR DEGREES
Overwrite
Echo Off
End
Molecule Cl-
Units BOHR
Cl 17.00 0.00000000 0.00000000 0.00000000 Type Cl
End
Read local pols for Cl-
Use ascii file Clminus_ref_wt4_L3_0f10.pol
Maximum rank 3
Limit rank to 3
Frequencies STATIC + 10
End
Molecule H2O
Units BOHR
O 8.00 0.00000000 0.00000000 0.00000000 Type O
H1 1.00 -1.45365196 0.00000000 -1.12168732 Type H
H2 1.00 1.45365196 0.00000000 -1.12168732 Type H
End
Read local pols for H2O
Use ascii file H2O_aTZ_ref_wt4_L3_0f10.pol
Maximum rank 3
Limit rank to 3
Limit rank to 1 for sites +++
H1 H2
Frequencies STATIC + 10
End
Write
File-prefix Clminus-H2O
CASIMIR file for Cl- +++
and H2O +++
with cutoff 0.0001
End
Finish
$ process < Clminus_H2O_casimir.prss > Clminus_H2O_WSM_L3iso.casimir
Title Cl- ... H2O
Frequencies 0.5 10
Skip 0
Print nonzero
Molecule Cl-
Site Cl type Cl
10 10
27.59015300 27.31627800 25.80430000 21.76549300 15.39099900
8.552807300 3.515195700 1.006707300 0.1707998800 0.6797508900E-02
11c 11c
27.59015300 27.31627800 25.80430000 21.76549300 15.39099900
8.552807300 3.515195700 1.006707300 0.1707998800 0.6797508900E-02
11s 11s
...
...
33s 33s
1094.410600 1088.040400 1050.364500 927.8203500 676.7228300
373.9877100 156.0034600 46.76507900 7.607165400 0.2612380000
End
Molecule H2O
Site O type O
10 10
6.483185500 6.461184100 6.331261300 5.904744500 4.952804400
3.465239500 1.847456100 0.6679903900 0.1282095700 0.5212038600E-02
11c 11c
...
...
Site H2 type H
10 10
1.427775400 1.419258000 1.370174900 1.221677900 0.9394025700
0.5798742900 0.2672248900 0.8300164000E-01 0.1333205500E-01 0.4515054600E-03
11c 11c
1.427775400 1.419258000 1.370174900 1.221677900 0.9394025700
0.5798742900 0.2672248900 0.8300164000E-01 0.1333205500E-01 0.4515054600E-03
11s 11s
1.427775400 1.419258000 1.370174900 1.221677900 0.9394025700
0.5798742900 0.2672248900 0.8300164000E-01 0.1333205500E-01 0.4515054600E-03
End
CGdir /home/alston/CamCASP/current/data/realcg
Dispersion 12 Cl- H2O
Finish
$ casimir < Clminus_H2O_WSM_L3iso.casimir > Clminus_H2O_WSM_L3iso.casimir.out
Dispersion coefficients for Cl- and H2O
Cl O C6 C7 C8 C9 C10 C11 C12
00 00 0 72.95569 0.0 1653.591 0.0 42806.65 0.0 858440.7
End
Cl H C6 C7 C8 C9 C10 C11 C12
00 00 0 14.14859 0.0 190.8377 0.0 2799.988
End
Basis Sets
Theory level
Energy evaluations with ORIENT