NumForce and symmetry

[onthefly]

Hi all,

When I used TM5.6, I was able to ensure the right symmetry to optimized excited states in the TD-DFT calculation of frequencies by NumForce. Now I'm using TM5.10, but I never succeeded in using C2v or Cs symmetry in NumForce TD-DFT excited states calculations.

NumForce is the TM5.10 tool to get  normal modes for ground state MP2,  CC2 ground/excited states or TD-DFT(HF) excited states. Actually TM5.10 manual says clearly that CC2 excited states must be treated C1 when using NumForce, and that with TD-DFT NumForce must be invoked with -n parameter whwre n is the n-th excited state of the 'a' irrep as molecule were C1 p.g.

Now the question. Apart from excited CC2 states, when is it possible to retain symmetry in NumForce calculations?

Thanks in advance.

onthefly

[christof.haettig]

I'm not sure what you mean with "ensure the right symmetry"...

NumForce does allways transform the input to C1 symmetry, because usually a large fraction of the displacements for the numerical differentiation break the point group symmetry. That's done independent of the method. The special problem for excited states is that the programs need to now the number of the excited state in C1.

Christof

[onthefly]

I'm thankful for the kind reply and I will try to be more specific. As an example consider trans glyoxal OHC-CHO, C2h symmetry. Well, I've been able to optimize its 1st excited state geometry (au irrep, ground is ag, of course) with TD-B3LYP / def2-SVP and def2-TZVPP basis sets, forcing glyoxal to C2h symmetry. Note that glyoxal has 3 unique atoms, so, in order to get normal modes with NumForce, only 18 SP/gradient calculations are needed. Now, when I run NumForce I always get the error attached below with SVP basis set; the same calculation indeed stops when I use TZVPP basis set and the only file TM creates in the Kraftwerk directory is the 'xm01.problem'. I had also a test with PBE0 functional obtaining the same results as B3LYP.

I checked all the 'coordinatenumber.coord' (.log, .dipole and .grad) files in the Kraftwerk directory and they are OK. When I don't use symmetry, (C1) with the same geometry (changes in the C1 optimized structure vs Cs structure are as small as 10^-6 a.u. )  I obtain all the normal modes with both functionals and basis sets, but now TM needs to execute 36 SP/gradient claculations and I have to assign normal mode symmetry 'by hand'.

Hope this will help to clarify my question.

Thanks

onthefly

@@@@@@@@@@@@@@@@@@@
forrtl: severe (64): input conversion error, unit 1, file /home/onthefly/glyoxal/S1/B3LYP/Cs/numforce/hessian
Image              PC        Routine            Line        Source
aoforce            08999794  Unknown               Unknown  Unknown
aoforce            0899733D  Unknown               Unknown  Unknown
aoforce            08967817  Unknown               Unknown  Unknown
aoforce            08933312  Unknown               Unknown  Unknown
aoforce            0894E51B  Unknown               Unknown  Unknown
aoforce            084A68D1  rdhss_.H                   59  rdhss.f
aoforce            080AB7D6  rdghd_.H                   36  rdghd.f
aoforce            08074C28  vibro_                    148  vibro.f
aoforce            080505C4  MAIN__.J                 1647  aoforce.f
aoforce            08048296  Unknown               Unknown  Unknown
aoforce            089A86A1  Unknown               Unknown  Unknown
aoforce            08048171  Unknown               Unknown  Unknown
complete results are in numforce/aoforce.out


                               Frequency Analysis

#############################################################
#                END OF NumForce                            #
#############################################################
 date: Wed Jul 30 16:28:56 CEST 2008

[acapobia]

I understand the problem about symmetry encountered by onthefly [who gave me his input files] because I experience it every time I run NumForce (TM5.10 parallel version), with both B3LYP (excited states) and CC2.

I don't know if it is actually a bug, but the 'xm01.problem' files comes from a wrong convergence of the SCF/SP related to the dscf calculation for the point of the 'xm01.coord' file.

I suppose it depends on the transformation of MOs from say C2h symmetry to C1 symmetry.  MOs are transformed the correct way, but the change of symmetry implies variations  in the 'metastase' 2e statistic file (for 2e integrals) also. However 'metastase' cannot be  updated if you have the  $statistics  off key in your control file.

Anyway, I have found a simple workaround.

Just run NumForce, then, when the 'KraftWerk' directory is created, just kill NumForce. Now go in the 'control' file of the /numforce dir and delete the lines

$statistics  off
$2e-ints_shell_statistics    file=metastase

then delete  the 'metastase' file.

Finally run NumForce again with the -c option; this will force dscf to make a statistic run each point (xm01.coord... zp03.coord and so on), and this run will generate a new metastase file, good for C1 symmetry, which (at my experience) leads to convergence for the dscf run and following egrad or ricc2 calculations.

Hope this will help.

Amedeo

[uwe]

Hi,

I am not sure where your problem really is, but there are three things that come to my mind:

• This is a point that you have already mentioned, but just for completeness:
  
  As Christof wrote: the problem is the NUMBER of the excited state you want to consider. NumForce does second numerical derivatives, so it transforms your system to C1 and shifts the atoms away a little bit from their initial positions.
  
  Now if you have C2h symmetry, you have Ag, Bg, Au, and Bu IRREPS, so you have chosen one of those excited states:
  
  irrep=au nexc=1
  
  i.e. the 1st excited state of the Au IRREP. You have to do a C1 single-point calculation to check which number that is in the numbering of C1's A IRREP.
  
  If you want to have a reduced number of steps, a possible work around would be to let C2h symmetry unchanged, do a short single-point check in C1 to get the number of the excited state you would like to get 2nd derivatives for, and then change the irrep=a nexc=<??> line for the C1 case.
  NumForce will transform the orbitals, but will not change the ricc2 keywords, so this should work.

• Please do not run the parallel version of Turbomole with NumForce!!!
  NumForce does a lot of independent single-point jobs, so if started with the -mfile option and the serial binaries, you will get ideal speed up - something you will not achieve when running all those jobs with the parallel version.
  In addition: NumForce is never being tested with the parallel binaries, so it is very likely that it will not work as one expects!

• Finally, there might be problems with the crossing of two near-lying excited states. Since one does a lot of single-point excited state calculations with slightly different coordinates, it can happen that two excited state change their position in the energy ranking.

Hope this helps a bit,

Uwe

[acapobia]

Uwe, thank you very much for the good tip of the -mfile option with the serial version as a workaround. I tested it and I obtained the same results as my workaround, however the procedure you suggest is considerably slower than running NumForce parallel version, at least on our cluster: moreover the problem of NumForce/parallel with symmetry at RI-CC2 level also affects  glyoxal (C2h) and e.g. acrolein (Cs)  ground states, where no flipping roots or crossing of states can occur.

Amedeo

[Hauke]

Hi,

I run into the same "problem" but am a bit surprised by the offered solution approach

let C2h symmetry unchanged, do a short single-point check in C1 to get the number of the excited state you would like to get 2nd derivatives for, and then change the irrep=a nexc=<??> line for the C1 case.
NumForce will transform the orbitals, but will not change the ricc2 keywords, so this should work.

Of course one has to change to "irrep=a" (actually, if one knows what one is doing it should be possible to use NumForce just for the preperation of the elongated structure files and then change the symmetry and irrep to the actual one of each elongated structure, submit them one by one and let NumForce just collect the results but I would consider this as advanced) but I would expect that one should specify the "number of the excited state you would like to get 2nd derivatives for"  by changing the state specified in the $ricc2 section "geoopt model=adc(2)    state=(a ??)". This section is most likely present from the structure optimization and the specified state might not exist in c1-symmtry (the usage of "state=(s1)" avoids editing this section but of course only works for the first excited state).

Maybe this was different in Turbomole 5 or am I missing some facts?