-
Notifications
You must be signed in to change notification settings - Fork 0
Tutorial
As an example I will estimate the static (hyper)polarizabilities of LiH at the CCSD/cc-pVTZ level of approximation using Gaussian'09 package. Assuming that python and numpy are properly installed and ff.py is in the search-path all I need is the file with Cartesian coordinates lih.xyz (in this case in atomic units).
2
Li 0.0 0.0 -1.50755
H 0.0 0.0 1.50755
Executing ff.py without any options returns the syntax:
This is a simple finite field helper script. By default it prepares a set of
input files for finite field calculations of each xyz structure provided as an
input, in which case an input template is required. If none is found a standard
one is prepared. Invoking the script with -c option allows to retrieve the
properties, provided that the calculations are completed, in which case paths
to data dir(s) are expected on input.
Usage: ff.py [options] xyz file(s) or data dir(s)
Options:
-h, --help show this help
-c, --calculate calculate properties for each data dir
-f, --base-field set base field ( >= 0.0001; defaults to 0.001 )
-s, --gaussian prepare gaussian inputs (expects gaussian.tmpl)
-g, --gamess prepare gamess inputs (expects gamess.tmpl)
-m, --molcas prepare molcas inputs (expects molcas.tmpl)
-u, --units set units; chose from: esu, si, asi
(si multiplied by electric permittivity of free space),
or au which is the default
--print-data if '-c' option is set, the script prints also the extracted
energies and dipole moments (usefull for testing purposes)
Fields selection:
--fields 25 25 fields sufficient for <g> (enabled by default)
13 13 fields sufficient for diagonal components
x 5 fields along x (y or z) axis for selected longitudinal components
GR[n,a,x] Use for Generalized Romberg: n - the number of iterations
a - the quotient (e.g. 2 for RR)
x - the axis (0,1,2 for x,y or z)
--print-fields print the chosen selection of fields
Machine epsilon is: 2.22044604925e-16 for float64 type
The following command should prepare a set of 25 input files (a minimum to get an isotropic average of the second hyperpolarizability), however, the script expects an input file template (in this case named gaussian.tmpl).
ff.py --gaussian lih.xyz
Since I did not provide any, a standard template should be punched out with the following info:
There's no gaussian template. I'm punching one - please check
OK. Now let us examine this default template file:
%chk=@chk
%mem=512mb
%nproc=2
#p hf/sto-3g scf(conver=11,novaracc,xqc) nosymm field=read
maxdisk=20gb iop(5/7=512) iop(3/27=20) iop(11/27=20)
integral(grid=ultrafine)
gaussian ffield
0 1
@data
@field
As you can see it is a bit paranoid, however, for reasonably accurate hyperpolarizabilities we do need highly accurate energies. The fields starting with '@' sign are going to be automatically updated. I need to change only the method, basis set and the default units to atomic:
%chk=@chk
%mem=512mb
%nproc=2
#p ccsd(conver=14,maxcyc=200)/cc-pVTZ scf(conver=11,novaracc,xqc) nosymm field=read
maxdisk=20gb iop(5/7=512) iop(3/27=20) iop(11/27=20) units=au
integral(grid=ultrafine)
gaussian ffield
0 1
@data
@field
This time, issuing the same command as before should prepare the required input files:
ff.py --gaussian lih.xyz
Now in the current directory I have the following [lih_0.0010_F00_.inp], [lih_0.0010_F01_.inp] and 23 similar files with different external field settings.
Now I'm ready to proceed with the calculations. Once all the calculations are finished I can use the ff.py again to perform the actual finite-field calculations. For that I need to run ff.py with -c option, providing the path to the directory with all the checkpoint files (*.chk) as an input. In this case I have all the files in the current directory ".".
ff.py --gaussian -c .
The script now attempts to execute the Gaussian utility formchk which I had to modify slightly in order to increase its default precision of output (see here for the particulars). Provided that the $g09root environment variable is set, the script should use the $g09root/g09/formchk binary, if not, it defaults to the /opt/gaussian/g09/formchk.
If everything goes smooth I should see the following messages for all the checkpoint files:
Using /opt/gaussian/g09/formchk to format checkpoints
Read checkpoint file ./lih_0.0010_F00_.chk
Write formatted file ./lih_0.0010_F00_.fchk
After all the formatted checkpoints are formed, the script extracts the necessary data and returns the following output for all the available energies (in the case of CCSD calculations these are: SCF, MP2, MP3, MP4(SDQ) and obviously CCSD). If the dipole moments are available ff.py calculates also the dipole-based properties.
Static electric dipole properties (finite field CCSD energy, f=0.0010)
Dipole moment [a.u.]
x y z <mu>
0.0000 -0.0000 -2.2926 2.2926
Polarizability [a.u.]
xx xy xz <alpha>
29.053 nan nan 27.386
yx yy yz
nan 29.053 nan
zx zy zz
nan nan 24.052
First Hyperpolarizability [a.u.]
xxx yxx zxx <beta_mu>
-0.00 0.00 442.84 -1043.42
xyy yyy zyy
-0.00 0.00 442.84
xzz yzz zzz
-0.00 -0.00 853.36
<beta_x> <beta_y> <beta_z> <beta>
-0.00 0.00 1043.42 347.81
Second Hyperpolarizability [a.u.]
xxxx xxyy xxzz <gamma>
47604.3 15859.2 30337.8 71619.8
yyxx yyyy yyzz
15859.2 47604.3 30338.0
zzxx zzyy zzzz
30337.8 30337.9 109820.4
The <beta_z> is the projection of the first hyperpolarizability tensor on the z-axis (and accordingly for <beta_x> and <beta_y>) given by eq. (5) in Ref. [1], whereas <beta_mu> is the similar projection on the direction of dipole moment. is simply the root mean square of the <beta_x>, <beta_y> and <beta_z>. and are the isotropic averages of the polarizability and second hyperpolarizability tensors assuming Kleinman's permutational symmetry (eqs. (4) and (6) of Ref. [1]).
Now I shall perform similar calculations using Gamess(US) package ... but to make them more interesting I shall attempt to estimate only the diagonal tensor elements using the more accurate Generalized-Romberg approach. Naturally I will use the same lih.xyz file and the calculations shall be performed using the same methodology i.e. CCSD/cc-pVTZ (actually this basis set is slightly different in both packages but the results should be virtually identical).
As before, I start by preparing a set of input files. Again, I need to prepare an input file template (in this case gamess.tmpl).
ff.py --gamess lih.xyz
Now a standard template should be punched with the following info:
There's no gamess template. I'm punching one - please check
Below is the standard Gamess(US) template with settings chosen to improve the accuracy of calculations.
$system mwords=10 memddi=0 parall=.t. $end
$contrl scftyp=rhf runtyp=ffield icharg=0 mult=1 units=angs
maxit=100 nprint=8 exetyp=run mplevl=0 ispher=-1
icut=20 itol=30 aimpac=.f. cctyp=none $end
!$eds mch(1)=0,0 mmul(1)=1,1 mnr(1)=1 ffeds=.t. $end
!$mp2 mp2prp=.t. code=ddi cutoff=1d-20 $end
$ffcalc offdia=.t. estep=@fstep onefld=@onefld
efield(1)=@field $end
!$efield evec(1)=@field sym=.f. $end
$ccinp iconv=14 maxcc=100 $end
$trans cuttrf=1d-15 $end
$scf dirscf=.t. fdiff=.f. diis=.t. soscf=.f.
conv=1d-11 swdiis=0.0001 $end
$basis gbasis=sto ngauss=3 $end
@data
I will change only the default model chemistry to CCSD/cc-pVTZ and the distance units to Bohr's. Moreover, I must switch off the parall option since the problem is simply to small to be parallelized.
$system mwords=10 memddi=0 parall=.f. $end
$contrl scftyp=rhf runtyp=ffield icharg=0 mult=1 units=bohr
maxit=100 nprint=8 exetyp=run mplevl=0 ispher=1
icut=20 itol=30 aimpac=.f. cctyp=ccsd $end
!$eds mch(1)=0,0 mmul(1)=1,1 mnr(1)=1 ffeds=.t. $end
!$mp2 mp2prp=.t. code=ddi cutoff=1d-20 $end
$ffcalc offdia=.t. estep=@fstep onefld=@onefld
efield(1)=@field $end
!$efield evec(1)=@field sym=.f. $end
$ccinp iconv=14 maxcc=100 $end
$trans cuttrf=1d-15 $end
$scf dirscf=.t. fdiff=.f. diis=.t. soscf=.f.
conv=1d-11 swdiis=0.0001 $end
$basis gbasis=cct $end
@data
This time, however, I will issue the following command in order to prepare the input files:
ff.py --gamess --fields GR[8,2,2] -f 0.0002 lih.xyz
The second option (--fields GR[8,2,2]) requests calculation of energies in the following uniform fields applied along z-axis: 0.0, +/- 0.0002, +/-0.0004, ..., +/-0.0256. That is +/-base field (0.0002) times the quotient (2) to the power of iteration no. (0,1,2,...,8). Now in the current directory I should find 17 files: lih_0.0002_F00_.inp, lih_0.0002_F01_.inp, ..., lih_0.0002_F16_.inp. Additionally, ff.py prepares lih_ffield_0.0002.inp for standard finite-field calculations as implemented in ''Gamess(US)'' package, but I will ignore this file.
Once all the calculations are successfully finished I can use the ff.py again to perform the analysis. Again, I need to run ff.py with -c option, providing the path to the directory with all the output files (*.log).
ff.py --gamess -c .
Please note that any additional file with *.log extension is likely to cause an error! This is because the script attempts to calculate all the possible properties based on the provided data. Since the names of the files are unimportant (only their content is) the results can be wrong also if the script finds two different files corresponding to the same external field. Therefore, it is recommended to collect only the required files in the working directory.
The output is now a bit more complicated since it contains not only the results of GR analysis but also the results estimated from the Kurtz formulas. For instance, for CCSD energies we find:
Static electric dipole properties (finite field CCSD energy, f=0.0002)
Dipole moment [a.u.]
x y z <mu>
nan nan -2.2926 nan
Polarizability [a.u.]
xx xy xz <alpha>
nan nan nan nan
yx yy yz
nan nan nan
zx zy zz
nan nan 24.052
First Hyperpolarizability [a.u.]
xxx yxx zxx <beta_mu>
nan nan nan nan
xyy yyy zyy
nan nan nan
xzz yzz zzz
nan nan 849.00
<beta_x> <beta_y> <beta_z> <beta>
nan nan nan nan
Second Hyperpolarizability [a.u.]
xxxx xxyy xxzz <gamma>
nan nan nan nan
yyxx yyyy yyzz
nan nan nan
zzxx zzyy zzzz
nan nan 109303.7
Since I have calculated only the energies for fields along z-axis, only the longitudinal components could have been calculated. However, this time I'm more interested in the results of Rutishauser-Romberg analysis (i.e. GR for the quotient equal to 2). These are printed in the form of the so-called Romberg's triangles for the 1st, 2nd, 3rd and 4th derivatives (for convenience multiplied by -1 in the case of energy based results). For example, in the case of 3rd derivative of CCSD energies (corresponding to the longitudinal first hyperpolarizability zzz) I have obtained:
GR scheme 3rd derivative (a=2.0) CCSD energy
-------------------------------------------------------------------------------
field / k 0 1 2 3 4 5 6
-------------------------------------------------------------------------------
0.0002 849.00 848.81 848.81 848.81 848.81 848.81 848.81
0.0004 849.55 848.82 848.82 848.82 848.82 848.82
0.0008 851.72 848.79 848.83 848.82 848.82
0.0016 860.52 848.23 848.94 848.78
0.0032 897.38 837.60 859.04
0.0064 1076.75 515.90
0.0128 2759.28
-------------------------------------------------------------------------------
P[p=2,k=3] = 848.82 +/- -0.00
------------------------------------------------------------------------------
The stability window in this case occurs for fields 0.0004-0.0008 after 2-5 iterations. The best estimate (corresponding to the smallest average deviations from the corresponding values in the previous iteration) is printed in the last row along with the rough estimate of numerical instability error. In this particular example one can see that the error in the case of standard finite-field calculations of the first hyperpolarizability (corresponding to the 0th iteration) with a typically assumed base field of 0.001 au would be about 3 au i.e. about 0.4%. The actual result of standard finite-field analysis (lih_ffield_0.001.log) with the base field of 0.001 au amounts to 853.36 au.
Since the results based on Gaussian'09 energies are slightly different, let us perform identical calculations in Gamess(US) package. First, I will make a new directory in which I will put the gamess template gamess.tmpl from previous calculations and lih.xyz file. Then I will prepare a set of required 25 files by issuing the following command in this directory:
ff.py --gamess lih.xyz
This is the default action corresponding to --fields 25 option. Also the base field of 0.001 au is the default setting. The results can be obtained after all the calculations had finished by executing ff.py with -c option:
ff.py --gamess -c .
The results below are virtually identical to those based on Gaussian'09 calculations with only small discrepancies in the case of second hyperpolarizabilities:
Static electric dipole properties (finite field CCSD energy, f=0.0010)
Dipole moment [a.u.]
x y z <mu>
-0.0000 -0.0000 -2.2926 2.2926
Polarizability [a.u.]
xx xy xz <alpha>
29.053 nan nan 27.386
yx yy yz
nan 29.053 nan
zx zy zz
nan nan 24.052
First Hyperpolarizability [a.u.]
xxx yxx zxx <beta_mu>
0.00 -0.00 442.84 -1043.42
xyy yyy zyy
0.00 0.00 442.84
xzz yzz zzz
0.00 0.00 853.36
<beta_x> <beta_y> <beta_z> <beta>
0.00 0.00 1043.42 347.81
Second Hyperpolarizability [a.u.]
xxxx xxyy xxzz <gamma>
47601.7 15857.5 30342.6 71622.0
yyxx yyyy yyzz
15857.5 47601.8 30342.7
zzxx zzyy zzzz
30342.6 30342.7 109820.9