SCREWER -- which "temperature" for normal mode elongation

[Jonas]

Dear community --

I frequently encounter saddle points of 1st (and eventually higher) order when carrying out molecular structure optimization calculations for some floppy metal ion--aqua complex systems. One possibility to come from saddle points to PES minima is to use SCREWER for elongation along the imaginary molecular normal modes found with frequency calculations. (I use, or have to use, NUMFORCE for that ...)
SCREWER asks for a "temperature", which appears to be a measure of normal mode elongation. However, I still feel unsure which "temperature" T to use, since a too low T won't get me far, and a too large T gets me too far from the PES minimum I feel I am quite close to ...

Is there any rule of thumb which T to use for given normal mode wave numbers?
Is there any one out there who has experience with the SCREWER tool? One does not find anything valuable on the WWW, I'm afraid ...

Any help or suggestions are welcome!
Thank you in advance,

Regards, Jonas

[Hauke]

It might be a bit late for Jonas but maybe my post is helpful for others who had the same question and find this thread.
By try and error I got to the following formula for the scaling factor that is used by screwer to shift the coordinates of the atoms.

(BoltzmannConstant*TempertureUsedInScrewer/ReducedMassOfModeOfInterest)^(1/2)/(2^(1/2)*Pi*SpeedOfLight*WavenumberOfModeOfInterest)

Here an example: For the 6th mode (the first one that is not 0) of CO₂  who has a wave number of 623.79 cm^-1  (I used just BP86 /def-SV(P) so don't wounder about the value) and a reduced mass (also written in output of aoforce, if you use NumForce look in numforce/aoforce.out) of 12.888 g/mol and 100 Kelvin used as temperature in screwer I get a factor of 0.0577696 Borrradius.

This factor is also written in the first line of the $newcoord section

Code: [Select]

$newcoord
# cartesian coordinates shifted along normal mode   6 by 0.05777
The atoms are then shifted by the product of this factor and the value that is given in the output of aoforce of this mode for each atom for x, y and z direction.

I tested this formula for H₂, D₂, CH₄ and CO₂ and it worked fine.

I also tried to understand this formula. If I didn't make a mistake it can be rewritten as:
(2/ForceConstantOfTheModeOfInterest*BoltzmannConstant*TempertureUsedInScrewer)^(1/2)

I don't know were this formula comes from (especially the factor of 2) but at least it makes sense (if the vibration has a huge force constant or if the temperature is low the shift is small).

But I agree that it is difficult to get a feeling which temperature is reasonable as the formula includes the reduced mass or force constant for which I have no "feeling" and scales only with the sqrt of the temperature. I would prefer to also have the choice to specify the scaling factor explicitly or use dimensionless displacement coordinates q...

[ole]

Thanks for your posting Hauke!

Your formula can be easily derived when starting from the (classical) harmonic oscillator potential E = 1/2 k x² resp. x = sqrt(2E/k). When you plug into that formula E=kT, k=4 pi² v² µ and v=cv^tilde, you arrive immediately at the term given. How did you find that by trial and error!? :-)

I "misused" the vibration program to do a potential energy scan along a vibrational mode (is there an easier way?). Now I wanted to plot the results as a function of the corresponding normal coordinate Q_i of the vibration, but got confused by the output of aoforce and the above calculated displacement factor.

In the aoforce output it says:
> each vibrational normal mode - given in terms of cartesian displacement vectors of all atoms - has been normalized to unity.
> to obtain mass-weigthed normal coordinates in a.u. divide the tabulated modes by sqrt(reduced mass * 1822.88853)
which implies that the displacement values given are not mass-weighted. But to get to mass-weighted normal coordinates, you have to multiply (and not divide) by sqrt(m_red) as in q = sqrt(m_red) x. Am I mistaken here? What is the dimension/unit of the displacements given by aoforce?

But then, it is easy to check that your formula for the shift is correct, a bit more formalized: x' = factor * deltaX * x. Therefore the shifting factor given by the vibration program should already be what you named "dimensionless displacement coordinate" Q_i. Am I missing a basic relationship here?

Actually I would expect all this information to be written in the manual, and not to be "reverse engineered" in the users forum :-( I'm not meaning to be rude!

It would be great if one of the developers would post here. Thank you!

Best, Ole

[christof.haettig]

Also the developers that work today on the package would need to get this info by 'reverse engineering'...

If someone has a good suggestion for a text which could be included in the manual, we will be happy to do so...

[JakubV]

Hello,

I would like to ask about this - we can approximate the PES in the direction of the imaginary freq. mode as

E(x) = E_0 - 1/2 kf x^2,

where E_0 is some constant, kf is a positive constant related to freq. as -i omega = sqrt(kf/m), where m is the reduced mass.
 
I wonder that the condition

E_0 - E(x) = kb T     (b for Boltzmann)

has two solutions for the new coordinates (here represented by the x).

x_1  =  x + deltax  (elongation)
x_2  =  x -  deltax (contraction)

(deltax = sqrt(2 kb T / kf), as was well derived and numerically tested above)

The screwer produce just the elongation one.

I am working on the SO4(2-).12H2O molecule (triplets of waters oriented around each sulfate oxygen, in the directions of lone electron pairs) - energy minimization by relax module and no symmetry considered.
I have got first time 4 rather small but imaginary frequencies, did screwing along each of the coordinates independently and the results had then just one or two im.freq., so I screwed them again, but even after iterating this procedure I still have one im.freq. left.
I have tried several temperature factors (from 50K to 350K), but I think the problem might that the mode is not as simple as for e.g. biatomic and the second solution may play a role.

 I can calculate it by hand, by the "screwer 2.0" should include the possibility for contraction calculation as well, I would suggest.

And, please would you suggest for my optimization problem?

I am using DFT-D3/B3LYP, def-SVP for all atoms (later I plan to use def-TZVPP for the sulfate part in order to better describe the least bonded electrons). My goal is to estimate (sulfate) complexation gibbs functions (and therefore stability constants) for several cations and lewis acids in general (with explicit hydration) - what do you think about this choice, please?


Best regards,
Jakub

[mpjohans]

I can calculate it by hand, by the "screwer 2.0" should include the possibility for contraction calculation as well, I would suggest.

How about entering a negative temperature?

Cheers,
    Mikael