TURBOMOLE Users Forum
TURBOMOLE Modules => Aoforce and Numforce => Topic started by: marand on February 17, 2021, 11:30:54 AM

I am trying to compute the lowest singlet excited state vibrations at the CC2/def2TZVPP level of theory for cisplanar conformer of bithiophene molecule.
The standard procedure I follow begins with geometry optimization for the excited state, then switching from C2v to C1 symmetry, having copied firs the gradient file
obtained during the geometry optimization into the submit directory.
The NumForce module (with level=cc2, central and ecnomic options) is then used and the calculations go smoothly till the hessian, aoforce output files are produced.
There is, however, one imaginary frequency that corresponds to the torsional motion of the rings along the central CC bond. This normally indicates a saddle point, which is in fact present in the ground state (in which goshconformers are the optimum structures). In the S1 state, however, increase coupling between the pi electrons of different rings makes the interring bond acquire a partly double character and rigidifies it, so the structure is expected to be actually planar. This has been reported in literature several times and my own experience also supports such a picture.
In order to check whether the imaginary frequency is some artifact of the numerical procedure or a genuine find, I have conducted DFT calculations (with B3LYP and camB3LYP functionals), which both returned allreal frequencies for the S1 excited states. Moreover, I have made a relaxed potential energy scan with respect to the interring rotation and it also shows that the planar conformation corresponds to the minimum on the PES. Finally, I used screwer to distort geometry in the direction of the imaginary frequency normal mode and then reoptimized the structure without any symmetry constraints (C1). The outcome was  again  a planar structure and the same imaginary frequency obtained for it in the NumForce calculations.
It, I think, confirms that the imaginary frequency is some numerical artifact. I tried to remove it by 1. increasing the convergence criteria for the orbitals, 2. changing the size of the geometry distortions in NumForce, 3. switching from ecnomic to cartesian displacements. None of those changed anything.
Please, advise me what can be the reason for such behaviour and how can I possibly get rid of the imaginary frequency. I need to interpret a helium temperature spectra and so I need vg quality of normal modes, including the lowfrequency ones.
Yours sincerely
Marcin Andrzejak

One more, but possibly important remark: I obtained the same results using the 7.0.1 and 7.5 versions of Turbomole.

Hi,
you mention that you run NumForce with "central" and "ecnomic" options. But aren't these in a sense conflicting options and if you give central before ecnomic, you may end up with the polyedr option. Could you confirm that you really do get "central" and not "polyedr" type displacements? I recommend running with central differences.
How large is the imaginary mode, by the way?
Best,
Antti

Well, the number of displacements shows that these are in fact central displacements, but I have already tried running the calculations without the ecnomic switch, and
the imaginary frequency is still there.
It is not large  i53 cm1, but it is nonetheless important for my studies to have it removed from the calculations.
Yours!
Marcin

Thanks for the clarification. Seems like a tough case. One question and one more comment:
1) You have used C_{2v} point group symmetry during the excited state optimization. Is there some particular reason why you remove symmetry before running NumForce? Is this because NumForce anyway breaks symmetry and then $excitations keyword used foir C_{2v} will not be valid? I'm just guessing here since I don't remember if I ever ran excited state frequency calculations with ricc2, but what if after the geometry optimization you just put in control file:
$ricc2
geoopt model=cc2 state=(s1)
$excitations
irrep=a nexc=1
# or possibly few more states in nexc, whatever works best for your system
And then run NumForce (for the C_{2v} symmetric structure). I'm just thinking if symmetry could help to reduce some numerical noise here.
2) You mentioned increasing convergence criteria for orbitals, but what keywords did you exactly use? Did you try something like
$scfconv 8
$denconv 1d8
or even tighter during the optimization? And then jobex with gcart 5?

Thanks for all the advice! However, I think I have tried those possibilities.
Geometry optimization was conducted in C2v symmetry, and only then the symmetry was artificially reduced to c1 so that the NumForce would run properly.
Coord and gradient was obtained with symmetry. The thight convergence for orbitals was required in exactly the way as shown, and the citeria were even tighter (10^10 for energy and density).
Geometry was optimized with the gcart 5 option.
The strangest thing is that the imaginary frequency appears even though the PES (along the imaginary frequency related normal mode) shows minimum, and geometry distorted along this coordinate returns to the planar one after reoptimization of the geometry. I must say, I am at a loss at to what to do now.
Best regards!
Marcin

Hi
I'm just bit confused about this one thing:
Geometry optimization was conducted in C2v symmetry, and only then the symmetry was artificially reduced to c1 so that the NumForce would run properly.
So did you manually reduce symmetry C2v > C1 before NumForce? Because this choice will significantly affect the number of displacements and possibly also have some numerical consequences (symmetry of the Hessian). What happens if you run NumForce directly for the C2vsymmetric structure?
Best,
Antti

I have never been able to compute excited state frequencies with symmetry. There has never been a complete hessian after all the displacement gradients were completed.
Ground state is ok, though.
What I normally do is I reduce symmetry to c1 (with define), set coordinates to cartesian, and finally change excitated state symmetry from e.g. b1 to a. So far I have never enountered problems of this kind, except of course close energetic proximity of some other state (which in c1 symmetry may interact with my target state), in which case all sorts of trouble may happen.
Can one do calculations of the excited state frequencies with symmetry, counting on the NumForce module to be smart enough to reduce the symmetry for the non totally symmetric displacements, changing the excited state representation accordingly? This would be marvellous, but I do not think it is ok with turbomole.
Yours!
Marcin

Hi,
About NumForce and symmetry: Earlier today I had a quick look and at least for a simple test system (C2v symmetry) this ran:
1. Using point group C2v, optimize the structure of the S1 state with jobex level cc2 (+other possible parameters)
2. Copy all files from the optimization directory to new directory.
3. Remove any C2vrelated ricc2 keyword from control (like $excitation_energies_CCS_____1^a1__ file=exstates)
4. Change $ricc2 data group to calculate gradients for S1 state in C1 point group:
$ricc2
geoopt model=cc2 state=(s1)
$excitations
irrep=a nexc=1
# or possibly few more states in nexc, whatever works best for your system
5. Run NumForce with
NumForce central c level cc2
It will generate displacements based on C2v. The gradients will be calculated for displaced structures (always with C1 symmetry and irrep a). Hessian will be obtained for C2vsymmetric.
So, in principle the process technically works, but I don't know if this is actually reasonable. H_{2}O runs fine, but I did get a large imaginary frequency. So, this is bit messy and I have doubts if this would actually help. Unfortunately, I don't have any other possible solutions in mind. Hopefully someone could also contribute!
Best,
Antti

Dear Antti, I believe I did just what you described. To be sure, I have just repeated the calculations, dnd the problematic imaginary frequency still remains.
I am quite positive that this is some kind of numerical artifact, but I do not know how can I get rid of it.
Perhaps it is possible to make Turbomole calculate gradients from 5 points instead of 3?
Marcin

Hi,
I'm not aware of any possibility to use 5 points instead of 3 in the numerical differentiation. The ruecker binary that generates the displacements only takes in $delta, $type (central/polyedr), $dimension (small, norm, big  this is a mystery to me!), and $freezeoption (nofreeze).
I agree that this sounds like a numerical artifact, the imaginary frequency is so small. From NumForce side, I only can think of two possibilities:
1) Using smaller or larger displacement (I guess you already tried these)
2) Using polyedr (if you so far used only central).
From ricc2 side, I don't know if increasing nexc could little bit affect the numerics. I guess you would see much larger imaginary frequency if this would be an issue of root flipping for some displacement.
I feel your pain, we are right now struggling with similar annoying small imaginary mode in phonon dispersion calculations for a solidstate material ;).
Best,
Antti

Hi,
just as a check: Distort the molecule along the vibrational mode with imaginary frequency and generate a series with different amplitudes (use TmoleX to do that or the vibration tool on the command line, run 'vibration help' first to see the options and then use the vibro option). Then, for each of the generated structures, run a singlepoint calculation and plot the energy (or excited state energy) vs the distortion. It will be very flat, but perhaps you can see whether it is a minimum or a maximum.
Regards,
Uwe

Thank you for a useful suggestion. I have done that and the PES scan says the molecule should be planar. although the surface is close to planarity in the vicinity of the supposed minimum.
However, I have also checked the augccpVTZ basis set and the addition of diffuse functions reduced the imaginary freqency for abou 50 to 7 cm1.
I guess I will have to live with that.
My very best regards!
Marcin