Isotope information ignored for vibrational/thermodynamic quantities

[TyBalduf]

Describe the bug
Isotopic information seems to have no effect on vibrational/thermodynamic quantities. From this commit I had thought isotope information would be incorporated.

To Reproduce
Steps to reproduce the behaviour: (I'm including ethanol as an example, but this has been the case for any molecule I have tried.

1. Run xtb ethanol.sdf --hess -v >> eth_regular.txt (vibrational calculation on regular ethanol)
2. Run xtb ethanol.sdf --hess -v --input temp.txt >> eth_iso.txt (vibrational calculation with all hydrogen mass set to 999.9 in the xcontrol) (isotope xcontrol keyword fails in a similar way)
3. Run xtb eth_iso.sdf --hess -v eth_iso_sdf.txt (vibrational calculation with a specific hydrogen mass set to 999.9 in the sdf

The three calculations give identical output, despite calculation (2) saying that the hydrogen was converted 999.9

Expected behaviour
Some difference between vibrational calculations with isotope information set. I don't know if xTB is supposed to be able to read this information from an sdf, but it seems like it should be able to from the xcontrol input.

Inputs: (Converted sdf to txt to upload)
eth_iso.sdf.txt
ethanol.sdf.txt
temp.txt

Outputs:
eth_iso.txt
eth_iso_sdf.txt
eth_regular.txt

[marcelmbn]

Confirmation/Reproduction of issue

Can confirm this, as I observed the exact similar problem without having seen this issue before.

Validation code

      call read_userdata(xcontrol, env, mol)
      write(*,*) "atmass in main l. 428", atmass
      ! ------------------------------------------------------------------------
      !> get some memory
      allocate (cn(mol%n), sat(mol%n), g(3, mol%n), source=0.0_wp)
      atmass = atomic_mass(mol%at) * autoamu ! from splitparam.f90
      write(*,*) "atmass in main l. 433", atmass

Validation output

mass of atom  1  changed to 400.00000000000000
 masses:   400.00000000000000        0.0000000000000000        0.0000000000000000        0.0000000000000000        0.0000000000000000        0.0000000000000000        0.0000000000000000        0.0000000000000000        0.0000000000000000        0.0000000000000000        0.0000000000000000        0.0000000000000000        0.0000000000000000        0.0000000000000000        0.0000000000000000     
mass of elements with Z=1 changed to 800.00000000000000
 masses:   400.00000000000000        0.0000000000000000        0.0000000000000000        0.0000000000000000        800.00000000000000        800.00000000000000        0.0000000000000000        0.0000000000000000        0.0000000000000000        800.00000000000000        0.0000000000000000        0.0000000000000000        800.00000000000000        800.00000000000000        800.00000000000000     
 atmass in main l. 428   400.00000000000000        0.0000000000000000        0.0000000000000000        0.0000000000000000        800.00000000000000        800.00000000000000        0.0000000000000000        0.0000000000000000        0.0000000000000000        800.00000000000000        0.0000000000000000        0.0000000000000000        800.00000000000000        800.00000000000000        800.00000000000000     
 atmass in main l. 433   14.006703209999998        12.010735900000000        12.010735900000000        12.010735900000000        1.0079407499999999        1.0079407499999999        15.999404920000000        15.999404920000000        14.006703209999998        1.0079407499999999        12.010735900000000        12.010735900000000        1.0079407499999999        1.0079407499999999        1.0079407499999999     

As you can see from the output, for this also quite simple organic molecule, I was able to temporarily change the atomic masses (in agreement with the output @TyBalduf mentioned, i.e. this is formally working). However, it's changed back due to this line: atmass = atomic_mass(mol%at) * autoamu ! from splitparam.f90. Thus, it's reverted to the physically correct ones. To me, it looks like just the ordering of the lines where physically correct values were initialized and user input is processed are interchanged.

I'll try to look a bit deeper into this.

In the later part of the numhess subroutine, consequently, only the physically correct values are used as we can see in this code/output:

Symmetrization and mass-weighting of the hessian

      ! symmetrize
      do i=1,n3
         do j=1,n3
            res%hess(j,i)=(h(i,j)+h(j,i))*0.5_wp
            if(abs(h(i,j)-h(j,i)).gt.1.d-2) then
               write(errStr,'(a,1x,i0,1x,i0,1x,a,1x,es14.6,1x,es14.6)') &
                  & 'Hessian element ',i,j,' is not symmetric:',h(i,j),h(j,i)
               call env%warning(trim(errStr), source)
            endif
         enddo
      enddo
      ! inverse mass array
      do ia = 1, mol%n
         write(*,*) "original atmass for index ",ia," is ",atmass(ia)
         do ic = 1, 3
            ii = (ia-1)*3+ic
            isqm(ii)=1.0_wp/sqrt(atmass(ia))
            amass(ii)=isqm(ii)/sqrt(amutoau)
        enddo
      enddo

           ------------------------------------------------- 
          |                Numerical Hessian                |
           ------------------------------------------------- 
step length          :   0.00500
SCC accuracy         :   0.30000
Hessian scale factor :   1.00000
frozen atoms in %    :   0.00000    0
RMS gradient         :   0.04898 !! INCOMPLETELY OPTIMIZED GEOMETRY !!
estimated CPU  time      0.01 min
estimated wall time      0.00 min
 original atmass for index            1  is    14.006703209999998     
 original atmass for index            2  is    12.010735900000000     
 original atmass for index            3  is    12.010735900000000     
 original atmass for index            4  is    12.010735900000000     
 original atmass for index            5  is    1.0079407499999999     
 original atmass for index            6  is    1.0079407499999999     
 original atmass for index            7  is    15.999404920000000     
 original atmass for index            8  is    15.999404920000000     
 original atmass for index            9  is    14.006703209999998     
 original atmass for index           10  is    1.0079407499999999     
 original atmass for index           11  is    12.010735900000000     
 original atmass for index           12  is    12.010735900000000     
 original atmass for index           13  is    1.0079407499999999     
 original atmass for index           14  is    1.0079407499999999     
 original atmass for index           15  is    1.0079407499999999     

[marcelmbn]

The issue was probably introduced in commit fa45f80 in the PR #930.

What was the reason for shifting the read_userdata within main.F90, @MtoLStoN, @demonic-daisy?

@Albkat, as you were responsible for the PR: We should look more closely into PRs which shift code within main.F90 where all the global variables are (unfortunately) modified on user-input.

[Albkat]

Thanks for sharing, we will definitely fix it in the upcoming release.

[MtoLStoN]

I have no idea anymore, this was more than a year ago. It was probably necessary to allow automated activation of cma constraints for the sandwich potential. As this should not interact with Isotope masses, shifting the assignment up should probably do the trick. Maybe it would also be a good idea to add a test case for this.

[marcelmbn]

I have no idea anymore, this was more than a year ago. It was probably necessary to allow automated activation of cma constraints for the sandwich potential. As this should not interact with Isotope masses, shifting the assignment up should probably do the trick. Maybe it would also be a good idea to add a test case for this.

Thanks for your reply and your insights (while you're busy somewhere else)! 💯

[MtoLStoN]

Sure thing. I think the userdata currently needs to be read before this line to allow the automatic cma setup for the sandwich potential.

xtb/src/prog/main.F90

Line 434 in c3cfd38

if (set%do_cma_trafo) then