TURBOMOLE Users Forum
TURBOMOLE Modules => Molecular Properties, Wavefunction Analysis, and Interfaces to Visualization Tools => Topic started by: prasanta13 on July 26, 2021, 02:21:23 PM
-
Hi there,
Is there any way to calculate the polarizability tensor for molecules?
Thanks in advance,
Prasanta Bandyopadhyay
-
Hello,
in the graphical user interface TmoleX you can select either static or dynamic polarizabilities in the Start-Job panel under the section "Spectra & Excited States". Then just run the job and after it finished you will directly see the tensor in the Results panel (the lower triangular part as the tensor is symmetric).
On the command line, for HF and DFT the keyword that has to be added to the input is $scfinstab (see chapter "Polarizabilities and Optical Rotations" in the documentation). For MP2 or CC2 the keyword $response is needed (see chapter "Ground State Second-order Properties with MP2 and CC2").
-
Thank You. It does help. ;D
-
I have successfully found the way to calculate the static polarizability with HF (from TmoleX) and generating control by hand.
However, I can not do the same for CCSD. TmoleX could not do that. I was trying to do that with define.
I could successfully implement ccsd within the ricc2 module, but within resp, there is only fop (first order property). I have also checked the manual "Ground State Second-Order Properties with MP2 and CC2". To compute the static orbital relaxed polarizabilities, I have to use,
$response
sop operators=(diplen,diplen) relaxed
Since, I could not find the the option from define, I added it manually. I also added the laplace=4 (without knowing what it does here).
Then I ran dscf module followed by ricc2 (which states that ccsd is in ccsdf12 , not in ricc2). Then I ran ccsdf12 which ended abnormally.
I am attaching the control and ccsdf12 output.
$title
$symmetry c1
$user-defined bonds file=coord
$coord file=coord
$optimize
internal off
redundant off
cartesian on
global off
basis off
$atoms
o 1 \
basis =o def2-TZVP
h 2-3 \
basis =h def2-TZVP
$basis file=basis
$scfmo file=mos
$scfiterlimit 30
$scfconv 7
$thize 0.10000000E-04
$thime 5
$scfdamp start=0.300 step=0.050 min=0.100
$scfdump
$scfintunit
unit=30 size=0 file=twoint
$scfdiis
$maxcor 500 MiB per_core
$scforbitalshift automatic=.1
$drvopt
cartesian on
basis off
global off
hessian on
dipole on
nuclear polarizability
$interconversion off
qconv=1.d-7
maxiter=25
$coordinateupdate
dqmax=0.3
interpolate on
statistics 5
$forceupdate
ahlrichs numgeo=0 mingeo=3 maxgeo=4 modus=<g|dq> dynamic fail=0.3
threig=0.005 reseig=0.005 thrbig=3.0 scale=1.00 damping=0.0
$forceinit on
diag=default
$energy file=energy
$grad file=gradient
$forceapprox file=forceapprox
$denconv 0.10000000E-06
$ricc2
ccsd
maxiter= 101
$laplace
conv=4
$response
sop operators=(diplen,diplen) relaxed
$rundimensions
natoms=3
nbf(CAO)=48
nbf(AO)=43
$closed shells
a 1-5 ( 2 )
$actual step ricc2
$orbital_max_rnorm 0.34522206215123E-05
$last SCF energy change = -75.668248
$charge from dscf
0.000 (not to be modified here)
$dipole from dscf
x 0.56305651579112 y 0.76490931398787 z 0.09993415747788 a.u.
| dipole | = 2.4274936608 debye
$end
The ccsdf12 outout:
data group $actual step is not empty
due to the abend of ricc2
this is only a subsequent error message
check reason for abend in the other output files ...
use the command 'actual -r' to get rid of that
========================
internal module stack:
------------------------
ccsdf12
turbo_setup
========================
CONTRL dead = actual step
ccsdf12 ended abnormally
What is the meaning of "data group $actual step is not empty
due to the abend of ricc2"?
What is/are the solutions?
-
This problem I have sorted out.
However, I tried to calculate static dipole polarizability with ccsd/def2-QZVPPD. I tried with the cc2 method also.
Below I attach the control file and output.
$title
$symmetry c1
$user-defined bonds file=coord
$coord file=coord
$optimize
internal off
redundant off
cartesian on
global off
basis off
$atoms
basis =def2-QZVPPD
cbas =def2-QZVPPD
$basis file=basis
$scfmo file=mos
$closed shells
a 1-5 ( 2 )
$scfiterlimit 30
$scfconv 7
$thize 0.10000000E-04
$thime 5
$scfdamp start=0.300 step=0.050 min=0.100
$scfdump
$scfintunit
unit=30 size=0 file=twoint
$scfdiis
$maxcor 500 MiB per_core
$scforbitalshift automatic=.1
$drvopt
cartesian on
basis off
global off
hessian on
dipole on
nuclear polarizability
$interconversion off
qconv=1.d-7
maxiter=25
$coordinateupdate
dqmax=0.3
interpolate on
statistics 5
$forceupdate
ahlrichs numgeo=0 mingeo=3 maxgeo=4 modus=<g|dq> dynamic fail=0.3
threig=0.005 reseig=0.005 thrbig=3.0 scale=1.00 damping=0.0
$forceinit on
diag=default
$energy file=energy
$grad file=gradient
$forceapprox file=forceapprox
$denconv 0.10000000E-06
$ricc2
cc2
$laplace
conv=4
$response
sop operators=(diplen,diplen) relaxed
$rundimensions
natoms=3
nbf(CAO)=158
nbf(AO)=132
$last step ccsdf12
$last SCF energy change = -.18016522E-09
$charge from dscf
0.000 (not to be modified here)
$dipole from dscf
x -0.00000000000002 y 0.67639454111482 z 0.38602176465619 a.u.
| dipole | = 1.9795155875 debye
$orbital_max_rnorm 0.24361361079503E-05
$cbas file=auxbasis
$end
Ran dscf and ccsdf12 module
I could not get the tensor, although in ccsdf12.out file. I think I needed to add more keywords in the control file which I can not figure out. From the previous reply and manual, I got this far.
This is my ccsdf12.out file
OpenMP run-time library returned nthreads = 48
ccsdf12 (beethoven) : TURBOMOLE rev. V7.5.0 compiled 13 Jun 2020 at 23:30:16
Copyright (C) 2020 TURBOMOLE GmbH, Karlsruhe
2022-03-08 12:04:41.380
************************************************************
* *
* C C S D F 1 2 P R O G R A M *
* *
* the quantum chemistry groups *
* at the universities in *
* Karlsruhe, Bochum, Bristol & Mainz *
* *
************************************************************
*-----------------------------------------------------------------------*
| program will use 48 thread(s) |
*-----------------------------------------------------------------------*
+--------------------------------------------------+
| Atomic coordinate, charge and isotop information |
+--------------------------------------------------+
atomic coordinates atom charge isotop
0.00000000 -0.10935034 -2.79640355 o 8.000 0
0.00000000 1.55511958 -3.50622571 h 1.000 0
0.00000000 0.15022504 -0.98141515 h 1.000 0
center of nuclear mass : 0.00000000 -0.00169876 -2.73456878
center of nuclear charge: 0.00000000 0.08305419 -2.68588692
+--------------------------------------------------+
| basis set information |
+--------------------------------------------------+
we will work with the 1s 3p 5d 7f 9g ... basis set
...i.e. with spherical basis functions...
type atoms prim cont basis
---------------------------------------------------------------------------
o 1 86 66 def2-QZVPPD [8s5p4d2f1g|16s9p4d2f1g]
h 2 36 33 def2-QZVPPD [4s4p2d1f|7s4p2d1f]
---------------------------------------------------------------------------
total: 3 158 132
---------------------------------------------------------------------------
total number of primitive shells : 46
total number of contracted shells : 42
total number of cartesian basis functions : 158
total number of SCF-basis functions : 132
symmetry group of the molecule : c1
the group has the following generators :
c1(z)
1 symmetry operations found
there are 1 real representations : a
=========================================================================
restricted closed shell calculation for the wavefunction models:
CC2 - Approximate CC Singles and Doubles
global parameters for ricc2 program:
hard restart (reuse of interm.) : disabled
soft restart (reuse of vectors) : disabled
threshold for vector function : 0.100000E-05
convergence threshold energy : 0.100000E-06
linear dependence threshold : 0.100000E-13
global print level : 1
maximum number of iterations : 30
maximum number DIIS vectors : 10
max. dim. of reduced space : 100
core memory limit (MB) : 24000
Scratch Directory :
=========================================================================
der. integral neglect threshold : 0.10E-07
integral neglect threshold : 0.21E-10
integral storage threshold THIZE : 0.10E-04
integral storage threshold THIME : 5
+------------------------------------------+
| Auxiliary basis set information |
+------------------------------------------+
assign orbital basis set names as defaults
read basis sets from /home/mainaks/turbomol/TURBOMOLE/cbasen/
(basis sets have been saved in file auxbasis)
we will work with the 1s 3p 5d 7f 9g ... basis set
...i.e. with spherical basis functions...
type atoms prim cont basis
---------------------------------------------------------------------------
o 1 157 157 def2-QZVPPD [10s9p8d6f3g1h|10s9p8d6f3g1h]
h 2 65 65 def2-QZVPPD [7s5p4d2f1g|7s5p4d2f1g]
---------------------------------------------------------------------------
total: 3 287 287
---------------------------------------------------------------------------
total number of primitive shells : 56
total number of contracted shells : 75
total number of cartesian basis functions : 373
total number of SCF-basis functions : 287
maximum number of shells which are related by symmetry : 1
The symmetry information takes 1 MB
*---------------------------------------------------------------------*
| simplified C1 algorithm will be applied |
*---------------------------------------------------------------------*
MOs are in ASCII format !
reading orbital data $scfmo from file mos
orbital characterization : scfconv=7
time elapsed for calculating density matrices : 0.316 sec
all orbitals will be active in the correlation treatment
Molecular Orbital Statistic:
============================
-----------------------------
orbitals in total:
-----------------------------
frozen occupied : 0
active occupied : 5
active virtual : 127
frozen virtual : 0
all together : 132
-----------------------------
time in riccmos cpu: 1.24 sec wall: 0.04 sec ratio: 28.0
total memory allocated for calculation of (Q|P)**(-1/2) : 1 MiB
calculation of (P|Q) ...
time in lp2sym cpu: 0.67 sec wall: 0.02 sec ratio: 28.0
calculation of the Cholesky decomposition of (P|Q)**(-1) ...
time in invmetri cpu: 0.18 sec wall: 0.01 sec ratio: 28.0
threshold for RMS(d[D]) in SCF was : 0.10E-06
integral neglect threshold : 0.21E-10
derivative integral neglect threshold : 0.10E-07
setting up bound for integral derivative estimation
increment for numerical differentiation : 0.00050000
=========================================================================
Energy of reference wave function is -76.0664790479200
Maximum orbital residual is 0.2436136107950E-05
Number of symmetry-nonredundant auxiliary basis functions: 287
Block lengths for integral files:
frozen occupied (BOI): 0 MiB
active occupied (BJI): 1 MiB
active virtual (BAI): 1 MiB
frozen virtual (BGI): 0 MiB
general (BTI): 1 MiB
=========================================================================
**************************************************************************
* *
* OPTIMIZATION OF THE GROUND STATE CLUSTER AMPLITUDES *
* *
**************************************************************************
start CC2 from scratch because restart not enabled
Iter. MP2 energy Norm(Omega) Norm(t1) Norm(t2) cpu wall
...........................................................................
1 -76.3760479185 0.0630498156 0.00000 0.15097 0.04 0.00
...........................................................................
time in total cpu: 2.72 sec wall: 0.06 sec ratio: 46.2
*****************************************************
* *
* RHF energy : -76.0664790479 *
* correlation energy : -0.3095688706 *
* *
* Final MP2 energy : -76.3760479185 *
* D1 diagnostic : 0.0156 *
* *
*****************************************************
Iter. CC2 energy Norm(Omega) Norm(t1) Norm(t2) cpu wall
...........................................................................
1 -76.3779715691 0.0094297715 0.02268 0.15284 0.06 0.00
2 -76.3785043180 0.0015036730 0.03256 0.15362 0.06 0.00
3 -76.3785218266 0.0002963720 0.03334 0.15364 0.06 0.00
4 -76.3785235950 0.0000438226 0.03355 0.15364 0.06 0.00
5 -76.3785238755 0.0000053247 0.03355 0.15364 0.06 0.00
6 -76.3785238955 0.0000006492 0.03355 0.15364 0.06 0.00
...........................................................................
CC equations converged in 6 iterations.
last energy change : 0.20E-07
convergence threshold : 0.10E-06
norm of the vector function: 0.65E-06
convergence threshold : 0.10E-05
time in total cpu: 22.09 sec wall: 0.47 sec ratio: 47.0
**********************************************************************
* *
* RHF energy : -76.0664790479 *
* correlation energy : -0.3120448476 *
* *
* Final CC2 energy : -76.3785238955 *
* *
* D1 diagnostic : 0.0251 *
* *
**********************************************************************
------------------------------------------------------------------------
total cpu-time : 37.94 seconds
total wall-time : 1.00 seconds
------------------------------------------------------------------------
**** ccsdf12 : all done ****
2022-03-08 12:04:42.330
Can anyone help me please?
-
Hello,
to calculate polarizabilities at CC2 level you need to add the right keywords to the input. Here is what the Turbomole documentations says:
[begin citation]
In addition to the standard input, second-order properties require that the data group for
the numerical Laplace transformation $laplace and that the sops option in the data group
$response is set. Frequency-dependent dipole polarizabilities with the CC2 model are
obtained with the input:
$ricc2
cc2
$laplace
conv=4
$response
sop operators=(diplen,diplen) freq=0.077d0
The frequency has to be given in atomic units.
Static orbital-relaxed polarizabilities are obtained with
$response
sop operators=(diplen,diplen) relaxed
[end citation]
So for static polarizability just add $ricc2 and $laplace as given above and $response as given in the second example.
Note that polarizabilites are implemented for CC2 and not for CCSD.
-
Thank you Uwe. Can you tell me if there is any upcoming plan to implement polarizabilities for CCSD and/or CCSD(T)?