ZPVC of Shielding not Independent of Molecule Orientation

[steska]

(Edit: In a reply I compiled data for a more basic example system, ethane)
Good day!

During the last few weeks I tried to do some research on the magnitudes and dependecies of zero-point vibrational corrections on NMR shieldings. I have, however, hit a wall when I realized the following:

A calculation's results are dependent upon the orientation of the input molecule.

The predicted shielding, corrected by the ZPV can be altered by up to 3ppm for a small molecule just by changing its input orientation.
As the effects I am looking for (environmental dependency of ZPVC) are supposedly in the area of 1ppm, research seems to be impossible at this point in time.

See the following example: A simple amide, called 1a (HCONH2) is used. The level of theory is HF/6-31G**. Dalton2015 is used. The input scripts are copypasted from the Dalton2015 manual and only slightly altered.

1. Geometric optimization of 1a.mol is performed using refgeo.dal
2. Minimum structure is extracted as 1a-ref.mol
3. Effective structures for 0K and 298K are calculated using effgeo.dal and extracted as 1a-eff000.mol and 1a-eff298.mol
4. 3 different rotations are generated for 1a-eff298.mol, rotating by 45° around each of the axes.
5. Shielding calculations are performed using the shield298.dal on the original and the three new orientations.
6. The vibrationally corrected shielding is calculated by taking the "Vibrationally corrected" values of Bxmx, Bymy and Bzmz from the "Nuclear magnetic shielding constants (ppm)"-subsection of the output and dividing by three.

The results (as an example for N):

Code: Select all

shield298_1a-eff298.out:     186.40225970
shield298_1a-eff298_x45.out: 186.00186449
shield298_1a-eff298_y45.out: 186.36358421
shield298_1a-eff298_z45.out: 188.29636215
Maximum difference:          2.3 ppm

All the while, the beginning of the output announces the structures to be the same (compare "Interatomic separations", starting line 370). The only thing different is the input orientation.

The same phenomenon appears when the rotation is applied to the initial structure or the reference structure.

What I played around with:

• All the threshholds and steplength
• .ACCURA (an undocumented five point derivation scheme that a colleague found while digging through the code
• different molecules
• sequential/parallel
• more stuff that I can't remember :p

Nothing had an effect on the described phenomenon.

To me it seems as though there are problems with the numerical approximations.
I also realized that the displacements (necessary for calculation of derivatives) are not performed along the normal modes but rather in x, y and z. Could this have an influnce?

Well, I'd be glad to get an answer, a reaction on this. Did I simply do something wrong? Underthunk something?

Thank your for reading! Have a nice day!
Best Regards,
Tim

[steska]

I have composed a simpler example with another basic system, ethane.

Assume the optimization has been carried out, as has been the generation of the effective geometries on HF/6-31g**-level.

The following effective structures of ethane for 298K have been generated by rotation:

Code: Select all

BASIS
6-31G**
Ethane rotated by 70 around y
Orientational depedancy of a pretty basic system
Atomtypes=2 Generators=0 
Charge=6.0 Atoms=2
C              -8.605376     0.000001    -5.635897
C              -5.874684    -0.000001    -4.642006
Charge=1.0 Atoms=6
H              -9.311148     1.917203    -5.892771
H              -8.743262    -0.958597    -7.452996
H              -9.879020    -0.958613    -4.332570
H              -5.736779     0.958614    -2.824917
H              -5.168922    -1.917205    -4.385125
H              -4.601054     0.958597    -5.945357

and

Code: Select all

BASIS
6-31G**
Ethane rotated by 70 around z
Orientational depedancy of a pretty basic system
Atomtypes=2 Generators=0 
Charge=6.0 Atoms=2
C               0.804704     2.210908   -10.013998
C               0.804705     2.210907    -7.108057
Charge=1.0 Atoms=6
H              -0.996879     2.866625   -10.765063
H               2.273366     3.443272   -10.765053
H               1.137637     0.322835   -10.765068
H              -0.663967     0.978565    -6.356987
H               2.606287     1.555181    -6.356999
H               0.471793     4.098992    -6.357008

Submitted through the copypasted code from the Dalton2015 manual:

Code: Select all

**DALTON INPUT
.WALK
*WALK
.VIBAVE
.DISPLACEMENT
0.05
.TEMPERATURES
1
298.0
**WAVE FUNCTIONS
.HF
*SCF INPUT
.THRESH
1.0D-10
**START
.SHIELD
*LINRES
.THRESH
1.0D-6
*RESPONS
.THRESH
1.0D-5
**EACH STEP
.SHIELD
*LINRES
.THRESH
1.0D-6
*RESPONS
.THRESH
1.0D-5
**END OF DALTON INPUT

The results are as follows:
The first part of the output (starting line 376) confirms that in fact, yes, we are looking at the same molecule

Code: Select all

  Bond distances (Angstrom):                                
  --------------------------                                
                                                            
                  atom 1     atom 2       distance          
                  ------     ------       --------          
  bond distance:  C          C            1.537758          
  bond distance:  H          C            1.089612          
  bond distance:  H          C            1.089613          
  bond distance:  H          C            1.089613          
  bond distance:  H          C            1.089613          
  bond distance:  H          C            1.089612          
  bond distance:  H          C            1.089613          
                                                            
                                                            
  Bond angles (degrees):                                    
  ----------------------                                    
                                                            
                  atom 1     atom 2     atom 3         angle
                  ------     ------     ------         -----
  bond angle:     C          C          H            111.393
  bond angle:     C          C          H            111.392
  bond angle:     C          C          H            111.393
  bond angle:     H          C          H            107.483
  bond angle:     H          C          H            107.483
  bond angle:     H          C          H            107.483
  bond angle:     C          C          H            111.393
  bond angle:     C          C          H            111.393
  bond angle:     C          C          H            111.392
  bond angle:     H          C          H            107.483
  bond angle:     H          C          H            107.483
  bond angle:     H          C          H            107.483

vs.

Code: Select all

  Bond distances (Angstrom):                                
  --------------------------                                
                                                            
                  atom 1     atom 2       distance          
                  ------     ------       --------          
  bond distance:  C          C            1.537758          
  bond distance:  H          C            1.089612          
  bond distance:  H          C            1.089613          
  bond distance:  H          C            1.089613          
  bond distance:  H          C            1.089613          
  bond distance:  H          C            1.089612          
  bond distance:  H          C            1.089613          
                                                            
                                                            
  Bond angles (degrees):                                    
  ----------------------                                    
                                                            
                  atom 1     atom 2     atom 3         angle
                  ------     ------     ------         -----
  bond angle:     C          C          H            111.393
  bond angle:     C          C          H            111.392
  bond angle:     C          C          H            111.393
  bond angle:     H          C          H            107.483
  bond angle:     H          C          H            107.483
  bond angle:     H          C          H            107.483
  bond angle:     C          C          H            111.393
  bond angle:     C          C          H            111.393
  bond angle:     C          C          H            111.392
  bond angle:     H          C          H            107.483
  bond angle:     H          C          H            107.483
  bond angle:     H          C          H            107.483

The vibrationally corrected shielding is calculated by adding up adding up the xx, yy and zz part of the result matrix and dividing by 3 (found around line 30970)
The results are:
shield298_00_00_70.out: C: 190.16499833 ; H: 30.75066417
shield298_00_70_00.out: C: 189.46655666 ; H: 30.76816624
Which, for Carbon, amounts to a difference of 0.5ppm

Any ideas?

[taylor]

I had intended to respond to your earlier posting but I was on travel. I do not understand your original statement that

The same phenomenon appears when the rotation is applied to the initial structure or the reference structure.

Are you implying that when you run a calculation of the chemical shifts at a fixed geometry (i.e., without any vibrational averaging) you still see a dependence on orientation? And/or that you see orientation dependence when you run with the temperature at 0 K?

I am puzzled too by your posted ethane geometries. The centre of mass seems to be considerably displaced from the origin of the coordinate system: how did this occur? That is, how did you arrive at these geometries in the first place? I am not asserting, by the way, that this should affect any calculations, but I wondered how it arose.

Finally, a point unrelated to your problem, if I may make it. 6-31G** is a reasonable choice if you are just running tests or trying to duplicate someone else's calculations. Otherwise, as a basis it is simply rubbish. Especially for properties. What rescues things, somewhat, for chemical shifts is that Dalton uses London orbitals and thus the basis set is perturbation-dependent for this particular property, which in turn helps accelerate convergence to the basis set limit in this case. So the results are merely not very good for this property, as opposed to completely worthless, which would be the case for many other properties. Just trying to be constructive here, by the way!

Best regards
Pete

[steska]

I had intended to respond to your earlier posting but I was on travel. I do not understand your original statement that

The same phenomenon appears when the rotation is applied to the initial structure or the reference structure.

What I meant is that when the rotation is performed before all the calculations (optimization, effective geometry, shielding) the problem still occurs. Just to make sure that Dalton isn't grabbing any old data without my knowledge and is then thrown off by the new orientation.

I am puzzled too by your posted ethane geometries. The centre of mass seems to be considerably displaced from the origin of the coordinate system: how did this occur? That is, how did you arrive at these geometries in the first place? I am not asserting, by the way, that this should affect any calculations, but I wondered how it arose.

You are right, I didn't really pay attention to that. It seems as if my very first geometry was already shifted (I just grabbed it out of chemcraft)

Code: Select all

BASIS
6-31G**
Ethane
Original, unoptimized structure
Atomtypes=2 Generators=0 Angstrom 
Charge=6.0 Atoms=2
C               1.245048     0.000000    -5.295295
C               1.245048     0.000000    -3.765307
Charge=1.0 Atoms=6
H               1.245048     1.019690    -5.694551
H               2.128126    -0.509845    -5.694551
H               0.361970    -0.509845    -5.694551
H               0.361970     0.509845    -3.366051
H               1.245048    -1.019690    -3.366051
H               2.128126     0.509845    -3.366051

You are also right of course, in that this should not really influence the calculations.

Finally, a point unrelated to your problem, if I may make it. 6-31G** is a reasonable choice if you are just running tests or trying to duplicate someone else's calculations. Otherwise, as a basis it is simply rubbish. Especially for properties. What rescues things, somewhat, for chemical shifts is that Dalton uses London orbitals and thus the basis set is perturbation-dependent for this particular property, which in turn helps accelerate convergence to the basis set limit in this case. So the results are merely not very good for this property, as opposed to completely worthless, which would be the case for many other properties. Just trying to be constructive here, by the way!

Just as HF is, as far as I am aware Yes, I only chose 6-31g** for testing and it would not have been the basis set used for actual property caclulation. I just ran some calculations using aug-cc-PVDZ, the difference in results depending on orientation is in the 0.5ppm area still.

Thank you for your interest and input!
Tim

[taylor]

If you make the effective temperature zero K is there still a dependence on orientation? And is it the case that the non-ZPVC results, i.e., the purely electronic contributions to the shift, are independent of the orientation? I am trying to see if we can narrow down which contributions are actually calculated as orientation-dependent.

By the way, in the original work on the ZPVC to properties we noted that in some cases the "vibrational frequencies" calculated at the effective geometry came out imaginary, particularly in the case of relatively low-frequency torsional modes. This should not as far as I can see cause orientation dependence, but it might be an issue.

Best regards
Pete

[taylor]

The answer may be simpler than all of this. If you compare your two posted outputs, it is true that the bond lengths and bond angles are the same for the two structures. However, this alone does not guarantee the structures are identical because the torsional angles might be different. And it seems they must be: if you compare the nuclear repulsion energies printed in the two outputs they are not the same. Therefore the two structures cannot be identical, differing only in orientation --- presumably however you have "rotated" the structure it must also have changed the torsional angle. There is then of course no reason why the chemical shifts should be the same in the two calculations.

Best regards
Pete
P.S. Early on we considered printing the torsional angles by identifying atoms that were appropriately connected, but the output rapidly grows and also tends to be highly redundant.

[taylor]

Can you please confirm that your "rotated" calculations are in fact not simply rotated but have changed the internal geometry of the molecule? We are approaching a new update release, and if there is actually a bug in the ZPV correction to properties we need urgently to know this and look to fixing it. It is my contention that your calculations are on two (slightly) different molecules and thus cannot be expected to yield the same results, in which case what you reported is not a bug (at least, not in Dalton!). It would help the Dalton community if we can resolve this, and expeditiously.

Best regards
Pete

[sauer]

Hi Pete,

I am sorry, but Tim cannot confirm anything right now, because he is already off on Xmas vacations back in Germany and I will go myself in a couple of days. But, when he comes back after New Year, we will look into your very useful comment. You are absolutely right, I should have asked him to check, whether the nuclear repulsion energy is the same.

Best wishes
Stephan

[steska]

Hey Pete and Stephan,

sorry, though I am already on my ways, I am still working. I will reply as soon as I have the time to come back to this, thanks for your efforts

[steska]

taking the results for two structures differing by 1.6 ppm

the torsional angles might be different.

This should be apparent from looking at the interatomic sperations:
1:

Code: Select all

           C           N           O           H           H           H       
           ------      ------      ------      ------      ------      ------  
C     :    0.000000                                                            
N     :    1.369949    0.000000                                                
O     :    1.192637    2.270597    0.000000                                    
H     :    1.102148    2.067756    2.010381    0.000000                        
H     :    1.992199    0.967950    2.471586    2.915292    0.000000            
H     :    1.999210    0.926055    3.078457    2.272693    1.675035    0.000000

2:

Code: Select all

           C           N           O           H           H           H       
           ------      ------      ------      ------      ------      ------  
C     :    0.000000                                                            
N     :    1.369949    0.000000                                                
O     :    1.192637    2.270597    0.000000                                    
H     :    1.102148    2.067756    2.010382    0.000000                        
H     :    1.992199    0.967950    2.471586    2.915292    0.000000            
H     :    1.999210    0.926055    3.078457    2.272694    1.675035    0.000000

which differ only in the last digit of one H--O and one H--H distance. A difference that small should not have an impact on calculation, should it?
As both molecules are recognized as being planar upon input, too, the torsional angles are the same as well.

if you compare the nuclear repulsion energies printed in the two outputs they are not the same.

1: @ Nuclear repulsion energy : 72.177600274522 Hartree
2: @ Nuclear repulsion energy : 72.177599669270 Hartree
difference: 6.05252e-7 Hartree
Is this difference of relevant magnitude? I compared the results for the 55 different rotations I generated and the maximum difference between two inputs in nuclear repulsion I could find was 1.3839858993947e-5 Hartree.

Best regards,
Tim

P.S.: Sorry for my late reply, I was quite sick :/

[taylor]

I am a bit confused now. The last discussions we were having were the simplified case --- ethane, where it is clear that the geometries are not identical: the nuclear repulsion energies differ by 2D-6 hartree. This is small, although larger than your new posting, where you have gone back, I guess, to your original molecule? Anyway, one's gut feeling would be such small differences should not be significant, but the only way to be sure is to eliminate any difference. If you are explicitly disabling symmetry, this should be easy --- run one calculation with the molecule in the xz plane and then swap the x and y coordinates so the molecule is now in the yz plane. This does not require any external program other than a text editor for the mol file. Does this give the same answer? If it does not, there is a bug. In principle all contributions should be rotationally invariant, at least for SCF (which you are using): the situation might be different for DFT. But if running exactly the same geometry (and there should be no problem in reproducing the nuclear repulsion energy exactly, at least to within machine precision) gives you the same results, then the problem lies elsewhere. And just confirming again, you are claiming that a single-point calculation at a particular geometry is not rotationally invariant even for the "static" (non-ZPV-corrected) property value?

Sorry to hear you've been unwell --- there's no good time to be ill bu Xmas/New Year must be one of the worst!

Best regards
Pete

[steska]

run one calculation with the molecule in the xz plane and then swap the x and y coordinates so the molecule is now in the yz plane.

As I did a lot of different rotations, this was of course one of them. If the molecule is rotated by 90° in any direction, the problem will not arise. The inbetweeners are where it gets messy.

And just confirming again, you are claiming that a single-point calculation at a particular geometry is not rotationally invariant even for the "static" (non-ZPV-corrected) property value?

I am not. The non-corrected values are the same.

Code: Select all

                                 At effective geometry:                         With ZPVC:
shield298_00_00_70.mol.out:           193.17817838                             190.16499833
shield298_00_70_00.mol.out:           193.17819383                             189.46655670

the nuclear repulsion energies differ by 2D-6 hartree. This is small, although larger than your new posting, where you have gone back, I guess, to your original molecule? Anyway, one's gut feeling would be such small differences should not be significant, but the only way to be sure is to eliminate any difference.

As with my original molecule, the structures in the ethane example differ in only two interatomic seperations by 1D-6Å each. The only way for me to make them even more similar (without using the 90° approach) is to turn up the significant digits. This should however not be necessary as we can see that the uncorrected shieldings are the same.

Could it be that the numerical differentiation is a) carried out along the 3 different axes and b) simply not accurate enough? No offense ment, just want to get input on this theory. This would also explain why the problem doesn't arise when rotating 90° and why the non-corrected shieldings don't show the problem.

Best wishes,
Tim

[taylor]

I am not convinced that we can make the statement

The only way for me to make them even more similar (without using the 90° approach) is to turn up the significant digits. This should however not be necessary as we can see that the uncorrected shieldings are the same.

Because there is as you note a difference between calculating the ZP corrections (using numerical integration) as compared to calculations at a static geometry. This is not an apples-to-apples comparison. So, without wanting to sound as though I keep playing the same tune over and over again: if you do increase the number of digits that specify the rotated geometry, does the discrepancy with the unrotated geometry decrease?

Best regards
Pete