Author Topic: Calculate the polarizability tensor for molecules  (Read 3408 times)

prasanta13

  • Full Member
  • ***
  • Posts: 36
  • Karma: +0/-0
Calculate the polarizability tensor for molecules
« 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

uwe

  • Global Moderator
  • Hero Member
  • *****
  • Posts: 558
  • Karma: +0/-0
Re: Calculate the polarizability tensor for molecules
« Reply #1 on: July 26, 2021, 03:30:50 PM »
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").

prasanta13

  • Full Member
  • ***
  • Posts: 36
  • Karma: +0/-0
Re: Calculate the polarizability tensor for molecules
« Reply #2 on: July 26, 2021, 04:18:15 PM »
Thank You. It does help. ;D

prasanta13

  • Full Member
  • ***
  • Posts: 36
  • Karma: +0/-0
Re: Calculate the polarizability tensor for molecules
« Reply #3 on: August 18, 2021, 09:28:30 AM »
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,
Code: [Select]
$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.

Code: [Select]
$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:
Code: [Select]
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?

prasanta13

  • Full Member
  • ***
  • Posts: 36
  • Karma: +0/-0
Re: Calculate the polarizability tensor for molecules
« Reply #4 on: March 08, 2022, 07:52:54 AM »
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.
Code: [Select]
$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
Code: [Select]
   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?

uwe

  • Global Moderator
  • Hero Member
  • *****
  • Posts: 558
  • Karma: +0/-0
Re: Calculate the polarizability tensor for molecules
« Reply #5 on: March 10, 2022, 09:55:19 AM »
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.




prasanta13

  • Full Member
  • ***
  • Posts: 36
  • Karma: +0/-0
Re: Calculate the polarizability tensor for molecules
« Reply #6 on: March 10, 2022, 01:16:19 PM »
Thank you Uwe. Can you tell me if there is any upcoming plan to implement polarizabilities for CCSD and/or CCSD(T)?