<> Navigation: * [[AJMPublic/camcasp/potentials|Potential Development using CamCASP]] * [[AJMPublic/camcasp|CamCASP]] * [[AJMPublic/orient|The Orient Program]] = 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 [[http://www-stone.ch.cam.ac.uk/documentation/gdma/|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. {{{#!wiki warning 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: 1. ''Distributed Multipole Analysis:  Stability for Large Basis Sets'', Stone, [[http://dx.doi.org/10.1021/ct050190+|J. Chem. Theory Comput., 1, 1128–1132 (2005)]]. This is the GDMA2 method for distributed multipoles. 1. ''Distributed polarizabilities obtained using a constrained density-fitting algorithm'', Misquitta and Stone, [[http://dx.doi.org/10.1063/1.2150828|J. Chem. Phys. 124, 024111 (2006)]]. These non-local polarizabilities form the basis for the WSM polarizabilities. 1. ''Accurate Induction Energies for Small Organic Molecules: 1. Theory'', Misquitta and Stone, [[http://dx.doi.org/10.1021/ct700104t|J. Chem. Theory Comput., 4, 7-18 (2008)]]. Describes the Williams-Stone-Misquitta (WSM) polarizability models 1. ''Accurate Induction Energies for Small Organic Molecules. 2. Development and Testing of Distributed Polarizability Models against SAPT(DFT) Energies'', Misquitta, Stone and Price, [[http://dx.doi.org/10.1021/ct700105f|J. Chem. Theory Comput., 4, 19-32 (2008)]]. Examples of the WSM polarizabilities. 1. ''Dispersion energies for small organic molecules: first row atoms'', Misquitta and Stone, [[http://dx.doi.org/10.1080/00268970802258617|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 [[AJMPublic/camcasp/polarization|Polarization]] and [[AJMPublic/camcasp/charge-transfer|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. {{{#!wiki important 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. {{{#!wiki important 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 == {{{#!wiki warning 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'''): [[attachment: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'''): [[attachment: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: 1. 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). 1. Run Cluster to generate the ORIENT input files for localization. 1. Run ORIENT to localize the non-local polarizabilities (using, in this example and by default, the Lillestolen-Wheatley localization scheme). 1. Run the Process code to write out the PFIT commands for the refinement stage. 1. Run PFIT to refine the localized polarizabilities. 1. Run Process to write out the Casimir commands. 1. 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'' {{{#!wiki important 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. }}} {{{#!wiki warning 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. {{{#!wiki important Important All results are in atomic units. }}} {{{#!wiki important 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 == {{{#!wiki warning 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'''): [[attachment:H2O-ISA.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'''): [[attachment: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 {{{ 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 }}} There may be other files too. The CamCASP output is in ''H2O.out''. Check and see what the ISA solution looks like. The ISA multipole moments are in ''H2O_ISA.mom'' and ''H2O_ISA-GRID.mom''. We recommend you use the latter as they are generally more accurate. {{{#!wiki important 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. }}} {{{#!wiki warning Warning What if the ISA does not converge? See [[https://app.ph.qmul.ac.uk/wiki/ajm:camcasp:isa#convergence_problems_with_the_isa|this page on possible solutions]] for help. }}} === Step (2) : ISA-Pol === Now we need to calculate the ISA-Pol non-local polarizabilities. This will be done as follows: * Create a directory for the ISA-Pol calculation. * Create a new Cluster file for the ISA-Pol calculation. * Execute Cluster to create a new CamCASP command file. * Link the MO file and the ISA solution file to the ISA-Pol directory. * Run CamCASP only using the 'runcamcasp' command (note, there is no ''.py'') 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 ''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. {{{#!wiki warning Warning The above commands assume you have a directory structure like this: H2O_aTZ/ /OUT/ /ISA-Pol/ 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 H2O/ /H2O_aTZ/ /OUT/ /ISA-Pol/ In this case the link commands would be: $ cd ISA-Pol $ ln -s ../H2O_aTZ/*movecs . $ ln -s ../H2O_aTZ/OUT/*ISA . }}} === Cluster file for ISA-Pol === The Cluster file for this calculation is very similar to the one use for the ISA calculation. There are some changes to the ''RUN-TYPE'' block that you will need to make: * The RUN-TYPE will be ''PROPERTIES'' * The AUX-BASIS will be the standard auxiliary basis with CARTESIAN GTOs. In the ISA calculation we used the ISA s-function set as part of this basis; here will will not do so. Why not? Because it does not seem to be needed. * Change the METHOD to ''isa-pol-from-isa-A-restart'' The complete Cluster file should look like this: [[attachment:H2O-ISApol.clt]] {{{ 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 }}} The axis file we will use is [[attachment: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 }}} {{{#!wiki warning 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. }}} === Running the ISA-Pol calculation === This time we do not need to run NWChem as we already have the molecular orbitals, so rather than use the ''runcamcasp.py'' - which cannot as yet skip the DFT part - we use an older Perl code which has the rather confusing name ''runcamcasp''. But as this script does not run Cluster, we need to do this first: {{{ $ cluster < H2O-ISApol.clt $ runcamcasp -q bg H2O_aTZ }}} {{{#!wiki warning 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 ''H2O_atoms.ISA''. Then the ISA-Pol calculation will start. When done you should have an ''OUT/'' directory containing the following files: {{{ 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 }}} The important files are: * H2O_ISA-GRID_f11_NL4_fmtB.pol : The ISA-Pol non-local polarizabilites in the new format. The older, format A, will soon be deprecated. * H2O_lim2.0_4.0_p2000_f11.p2p : point-to-point polarizabilities. * H2O_aTZ.out (or H2O_aTZ_camcasp.out): CamCASP output file. Look at this to see if the results make sense. The total polarizabiities and some of the non-local polarizabilities are printed in this file. {{{#!wiki warning 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: [[attachment:H2O total polarizabilities]] {{{ 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 }}} These were calculated with a larger basis, but the aug-cc-pVTZ basis should result in something like this. === Localization of the ISA-Pol polarizabilities === The non-local polarizabilities can be localized using the ''localize.py'' code. This is a multi-step process which has been described above. Before performing the localization we need to make some decisions: * The maximum rank of the localized polarzabilities. * Is this to be an isotropic or anisotropic model? At present we do not have the ability to generate hybrid models, but this will be done soon (it can be done with a manual editing of the intermediate files created in the localization process). * The local axis frame for the anisotropic model. We have already defined the axis frame for H2O (see above). * Some settings used in the refinement stage. Let's have a look at both isotropic and anisotropic models. === Isotropic local models === Here we will create a rank 3 isotropic local model for water. Make sure you are in the ''ISA-Pol/'' directory and issue the following command: {{{ localize.py --limit 3 --wsmlimit 3 --hlimit 3 H2O_aTZ }}} 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: [[attachment:loc-iso-command.bash]] {{{ localize.py --limit 3 --wsmlimit 3 --hlimit 3 \ --isotropic --subdir L3iso-wt3-coeff1e-3 \ --loc LW --weight 3 --weightcoeff 1e-3 \ H2O_aTZ }}} {{{#!wiki important 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: $ chmod +x loc-iso-command.bash Now you can execute the commands in this file using: $ ./loc-iso-command.bash }}} Here are the options used in the above: * --limit 3 : perform the initial localization to rank 3. * --wsmlimit 3 : Set the rank of the refined polarizabilities to 3. This applies to all but hydrogen atoms, for these use * --hlimit 3 : to set the rank of the polarizabilities on hyrdogen atoms also to 3. * --isotropic : force isotropic models * --subdir L3iso-wt3-coeff1e-3 : do the calculations in sub directory ''L3iso-wt3-coeff1e-3/'' * --loc LW : use the Lillestolen-Wheatley localization for the initial step. * --weight 3 : use weight scheme 3 for the refinement step * --weightcoeff 1e-3 : and associated with this weight is a strength, set this to 0.001 in arbitrary units. This is fairly strong and it keeps the polarizabilities from wandering too far from the initial values. * H2O : the file prefix. When executed you will get the following output: {{{ $ ./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. }}} {{{#!wiki warning 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: [[attachment: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 }}} One way of assessing if this is a good model is to have a look at the output files of the form ''H2O_ref_*_xxx.out''. The one with ''xxx=000'' contains the output of the refinement of the static polarizabilities. It's quite a large file, but the important parts are: {{{ ... ... 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 }}} Look for the following: * The residuals: they should be small. 6% is OK for an isotropic model. * The total polarizability: This should be close to what was calculated by CamCASP. === Anisotropic local models === The only difference with anisotropic models is that you need to have the axis file present in the work directory (''ISA-Pol/''), and that the refinement takes a lot longer as there are more parameters to optimise. {{{#!wiki warning Warning Make sure you have copied or linked the ''H2O.axes'' file to the ''ISA-Pol/'' directory! }}} Here is the localization command: [[attachment:loc-aniso-command.bash]] {{{ localize.py --limit 3 --wsmlimit 3 --hlimit 3 \ --axes H2O.axes --subdir L3-wt3-coeff1e-3 \ --loc LW --weight 3 --weightcoeff 1e-3 \ H2O }}} Execute this as before... {{{ $ ./loc-aniso-command.bash }}} And the output will be in directory ''L3-wt3-coeff1e-3/''. The potential file now looks like this: [[attachment:H2O_ref_wt3_L3_C12.pot]] {{{ ! 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 ... ... }}} Once again we should assess the localization by looking at the file ''H2O_ref_wt3_L3_000.out'': [[attachment:H2O_ref_wt3_L3_000.out]] {{{ ... ... 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 }}} 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). {{{#!wiki warning 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. * Your axis file may not be in the work directory. }}} == Hetero-dimers : Dispersion models == 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). * The Process command files that create the Casimir input files: These are //Clminus_casimir.prss// and //H2O_aTZ_casimir.prss//. 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: {{{ 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 }}} 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: {{{ $ process < Clminus_H2O_casimir.prss > Clminus_H2O_WSM_L3iso.casimir }}} This file looks like: {{{ 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 }}} {{{#!wiki important 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: {{{ $ casimir < Clminus_H2O_WSM_L3iso.casimir > Clminus_H2O_WSM_L3iso.casimir.out }}} The output contains the dispersion model. It will be at the end of this file. For this example we get: {{{ 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 }}} {{{#!wiki important Important If anisotropic polarizabilities are used, the resulting dispersion model will be in the axis system used for those polarizabilities! }}} == Basis Sets == 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. == Theory level == 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. == Energy evaluations with ORIENT == 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 [[ajm/orient/energy-calculations|Orient: Energy Calculations]] section for more details.