Problem with DFTB+ when using i-PI

[Nicolas Künzel]

Hello,

right now I am trying to set up a calculation with DFTB+ and i-PI but I experience some problems.
I take an input-file (dftb_in.hsd) with a .gen-file (bulkpair.gen) which works well when one uses the built in Velocity-Verlet driver of DFTB+. But when I change the driver to i-PI (as shown in the dftb_in.hsd file) then the parser starts and works well until some point where I get an error. The whole output is:

***  Parsing and initializing

Parser version: 5

Interpreting input file 'dftb_in.hsd'
--------------------------------------------------------------------------------
***  Converting input from version  4 to version  5 ...
WARNING!
-> Keyword renamed to 'QR'.
Path: dftb_in/Hamiltonian/DFTB/Eigensolver/Standard
Line: 28-28 (File: dftb_in.hsd)

WARNING!
-> Adding legacy step size for finite difference differentiation
Path: dftb_in/Hamiltonian/DFTB/Differentiation

***  Done.

Reading SK-files:
  /cluster/apps/dft/Slater-Koster-Lib/3ob-1-1/C-C.skf
  /cluster/apps/dft/Slater-Koster-Lib/3ob-1-1/C-N.skf
  /cluster/apps/dft/Slater-Koster-Lib/3ob-1-1/C-H.skf
  /cluster/apps/dft/Slater-Koster-Lib/3ob-1-1/C-O.skf
  /cluster/apps/dft/Slater-Koster-Lib/3ob-1-1/N-C.skf
  /cluster/apps/dft/Slater-Koster-Lib/3ob-1-1/N-N.skf
  /cluster/apps/dft/Slater-Koster-Lib/3ob-1-1/N-H.skf
  /cluster/apps/dft/Slater-Koster-Lib/3ob-1-1/N-O.skf
  /cluster/apps/dft/Slater-Koster-Lib/3ob-1-1/H-C.skf
  /cluster/apps/dft/Slater-Koster-Lib/3ob-1-1/H-N.skf
  /cluster/apps/dft/Slater-Koster-Lib/3ob-1-1/H-H.skf
  /cluster/apps/dft/Slater-Koster-Lib/3ob-1-1/H-O.skf
  /cluster/apps/dft/Slater-Koster-Lib/3ob-1-1/O-C.skf
  /cluster/apps/dft/Slater-Koster-Lib/3ob-1-1/O-N.skf
  /cluster/apps/dft/Slater-Koster-Lib/3ob-1-1/O-H.skf
  /cluster/apps/dft/Slater-Koster-Lib/3ob-1-1/O-O.skf
Done.


Processed input in HSD format written to 'dftb_pin.hsd'

Starting initialization...
--------------------------------------------------------------------------------
 Initialising for socket communication to host /tmp/ipi_dftb_equilibration
Connecting to :/tmp/ipi_dftb_equilibration:
Mode:                        Socket controlled calculation
Self consistent charges:     Yes
SCC-tolerance:                 0.100000E-08
Max. scc iterations:                   1000
Ewald alpha parameter:         0.000000E+00
Spin polarisation:           No
Nr. of up electrons:            34.000000
Nr. of down electrons:          34.000000
Periodic boundaries:         No
Diagonalizer:                Standard
Mixer:                       Anderson mixer
Mixing parameter:                  0.050000
Maximal SCC-cycles:                    1000
Nr. of chrg. vectors to mix:              8
Force evaluation method:     Traditional                                                                                                                                                                                            
Electronic temperature:        0.100000E-07
Initial charges:             Set automatically (system chrg:   0.000E+00)
Included shells:             C:  s, p
                             N:  s, p
                             H:  s
                             O:  s, p
Using Lennard-Jones dispersion corrections
Full 3rd order correction    Yes
Damped SCC                   Yes
Damped species(s):           H                                                
Extra options:
                             Mulliken analysis
Force type                   original
 

Intel MKL ERROR: Parameter 4 was incorrect on entry to DGETRF.
ERROR!
-> Failure in LU factorisation dgetrf, illegal argument at position         -4


As I found on Google that is a lapack/blas error.
So maybe you already now that error and can tell me what I can do there in order to get it working?

Cheers
Nicolas

[Michele Ceriotti]

Does the DFTB+ exampled work (e.g. the Zundel)? Generally this kind of errors might come from crazy atomic configurations being passed on, which might be due to the i-PI input being provided (or interpreted) in the wrong units.

Michele

[Nicolas Künzel]

Actually the only DFTB+ examples I could find in i-PI are outdated. The one in the master branch doesn't even include an input file for DFTB+ and the one in the dftb-rem-bias branch have input.xml files which are not working with the recent version of i-PI.

[Michele Ceriotti]

We forgot adding the DFTB+ inputs. Pull from master, now they're in.

[Nicolas Künzel]

Ok. So I got the files and I think some things have to be changed there.

First: The .hsd file should be named dftb_in.hsd and not dftb_pin.hsd (pin is the already parsed file which doesn't work as direct dftb+ input, dftb+ is doing that by itself).
Second: In the .hsd file the driver should be Socket and not i-pi. So there is smth wrong.

And after changing all of that I get the following error:

***  Parsing and initializing

Parser version: 5

Interpreting input file 'dftb_in.hsd'
--------------------------------------------------------------------------------
***  Converting input from version  4 to version  5 ...
WARNING!

-> Keyword moved to Analysis block.
Path: dftb_in/Options/AtomResolvedEnergies
Line: 85-85 (File: dftb_in.hsd)

WARNING!
-> Keyword moved to Analysis block.
Path: dftb_in/Options/WriteEigenvectors
Line: 82-82 (File: dftb_in.hsd)

WARNING!
-> Keyword moved to Analysis block.
Path: dftb_in/Options/WriteBandOut
Line: 81-81 (File: dftb_in.hsd)

WARNING!
-> Adding legacy step size for finite difference differentiation
Path: dftb_in/Hamiltonian/DFTB/Differentiation

***  Done.

Reading SK-files:

  ../dftb-param/O-O.skf
  ../dftb-param/O-H.skf
  ../dftb-param/H-O.skf
  ../dftb-param/H-H.skf
Done.
WARNING!
-> The following 1 node(s) have been ignored by the parser:
(1)
Path: dftb_in/Hamiltonian/DFTB/OldRepulsiveSum
Line: 59-59 (File: dftb_in.hsd)

ERROR!
-> Code halting due to the presence of errors in dftb_in file.

[Nicolas Künzel]

Ok when removing the OldRepulsiveSum entry from the .hsd file it works well.

[Nicolas Künzel]

So the connection between i-PI and DFTB+ works well.

And when I use my own input just with DFTB+ and do the MD in there it works as well.
But when I take my input, use it with i-PI and DFTB+ I get this weird error.
Could it be a problem that i-PI is not working with the "cluster"-format (non-periodic) of DFTB+ or smth?
I really hoped I could just take the input as it is from another guy from our group and try to get the same results he got with another method.

[Michele Ceriotti]

The problem (as you can see from the header) is that these inputs were made for the very first patched version of dftb+, which was 4.0
So the inputs need to be updated to v5 -- the examples run perfectly on my laptop using the patched v 4.0.
So this actually is more of a dev issue: can you please open a bugreport on github, and perhaps help out making the example DFTB+ work with v5.0?
Then we can merge. I suspect your own problems, instead, come from an issue with units or with the interpretation of the input (cell, etc.).
Look at the first xyz file that should be printed reflecting the starting configuration as read by i-PI. Does that make sense?
Best
Michele

[Michele Ceriotti]

Sorry I had missed your updated messages. So one possibile problem is that i-PI by defaults wraps atoms in the unit cell before passing them to the driver. This comes because some codes (e.g. LAMMPS) expect internal coordinates to be close enough to being within the unit cell.  You can disable this by setting the attribute pbc="False" in the <ffsocket> definition.

[Nicolas Künzel]

So I did both, print out the trajectory every time step so I get the starting configuration. It looks the same as the input.
The xyz-file I put in was:

   27
#CELL{H}: 71.230 0.0 0.0 0.0 71.230 0.0 0.0 0.0 71.230 positions{angstrom}
    N      7.92759019      1.88677257      6.66054458
    C      7.58554909      0.67905532      5.71618874
    C      9.22458898      2.51324801      6.32301501
    C      6.80492718      3.07210038      6.67452946
    C      8.10866699     -0.63310037      6.24253551
    C      9.84720867      3.11567471      7.57739865
    C      6.28801226      3.40495043      5.30758164
    H      7.83577028      1.52708928      7.79657812
    H      6.45991992      0.60924005      5.69495240
    H      7.96318445      0.93669442      4.72556261
    H      9.85403940      1.76773747      5.80892846
    H      8.97597135      3.33986945      5.67383631
    H      7.23717042      3.95170768      7.16546316
    H      5.95694235      2.83673752      7.31731158
    H      8.03092080     -1.40077791      5.50319067
    H      7.52987039     -0.99028972      7.04382312
    H      9.13561828     -0.61893109      6.52444501
    H      9.90842061      2.37924447      8.38793062
    H     10.88280016      3.47580700      7.27279664
    H      9.22928755      3.96639021      7.87292968
    H      5.67856003      2.61067064      4.85064503
    H      5.63769852      4.26540913      5.43794499
    H      7.12650903      3.64634250      4.64773905
    O      9.51585820      0.39085824      9.29723553
    O      8.24310437      0.63094721     11.13422483
    O      7.49968532      1.26570472      9.15232222
    N      8.43941916      0.67477877      9.92092336


and the printed out one is:

27
# CELL(abcABC):   71.23000    71.23000    71.23000    90.00000    90.00000    90.00000  cell{atomic_unit}  Traj: positions{angstrom} Step:           0  Bead:       0
       N  7.92759e+00  1.88677e+00  6.66054e+00
       C  7.58555e+00  6.79055e-01  5.71619e+00
       C  9.22459e+00  2.51325e+00  6.32302e+00
       C  6.80493e+00  3.07210e+00  6.67453e+00
       C  8.10867e+00 -6.33100e-01  6.24254e+00
       C  9.84721e+00  3.11567e+00  7.57740e+00
       C  6.28801e+00  3.40495e+00  5.30758e+00
       H  7.83577e+00  1.52709e+00  7.79658e+00
       H  6.45992e+00  6.09240e-01  5.69495e+00
       H  7.96318e+00  9.36694e-01  4.72556e+00
       H  9.85404e+00  1.76774e+00  5.80893e+00
       H  8.97597e+00  3.33987e+00  5.67384e+00
       H  7.23717e+00  3.95171e+00  7.16546e+00
       H  5.95694e+00  2.83674e+00  7.31731e+00
       H  8.03092e+00 -1.40078e+00  5.50319e+00
       H  7.52987e+00 -9.90290e-01  7.04382e+00
       H  9.13562e+00 -6.18931e-01  6.52445e+00
       H  9.90842e+00  2.37924e+00  8.38793e+00
       H  1.08828e+01  3.47581e+00  7.27280e+00
       H  9.22929e+00  3.96639e+00  7.87293e+00
       H  5.67856e+00  2.61067e+00  4.85065e+00
       H  5.63770e+00  4.26541e+00  5.43794e+00
       H  7.12651e+00  3.64634e+00  4.64774e+00
       O  9.51586e+00  3.90858e-01  9.29724e+00
       O  8.24310e+00  6.30947e-01  1.11342e+01
       O  7.49969e+00  1.26570e+00  9.15232e+00
       N  8.43942e+00  6.74779e-01  9.92092e+00

Do you see smth weird?

I also added the pbs="False" in <ffsocket>, that didn't change anything either. I am a bit out of ideas right now.

[Nicolas Künzel]

So I found the error. When I change the input for DFTB+ in the gen format from cluster (C) to supercell (S) then it works. Maybe i-PI sends cell data and if that is not specified in DFTB+ there are errors?

[Michele Ceriotti]

This sounds like a problem in the DFTB+ driver implementation - the (minimalistic) i-PI protocol always passes cell size so perhaps these should just be ignored if the geometry is a cluster one. I added Ben to this list so perhaps he can comment and if we agree this is a driver-side problem, correct it. 

[Ben Hourahine]

Hi,

the master branch on the DFTB+ github (https://github.com/dftbplus/dftbplus) can cope with cluster geometries in the input while using a socket interface (lattice vectors passed by i-PI are ignored in this case, so beware if i-PI is assuming it can barostat or whatever). Unfortunately the last DFTB+ release (17.1) assumed that since the i-pi protocol sends lattice information the geometry would always be periodic (and didn't trap this correctly).

If the calculation is non-self-consistent, just setting a large unit cell would be equivalent to cluster geometries. In charge self-consistent mode (SCC=Yes) you can use the same strategy, but there are two things to be aware of. Firstly the Ewald summation parameters in the main code can become inefficient in very large (or irregular) shaped periodic systems (the next 17.2 release will do something about this). Secondly, the electrostatics for a periodic system cause there to be some issues related to the total dipole moment of the unit cell being forced to be 0, however normally you don't see these artefacts in real calculations.

Regards

Ben

[Nicolas Künzel]

Thanks a lot for the answers!

I guess I will do it as written by Ben then.