strange curve of the binding energy of Na-Cl

[mejdeddine mokhtar]

Dear CP2K developers and users,

I've tried to compute the binding energy of Na _ Cl (singlet state so the Multiplicity should be 1)but I got a strange curve at large distances. I don't have a lot o experience in using cp2k so I attached the script that I've used to calculate the potential energy at varying distance. Could someone check if I did something wrong? any help would be appreciated. 

&GLOBAL

    PROJECT NaCl

    RUN_TYPE ENERGY             

&END GLOBAL

&FORCE_EVAL

  METHOD Quickstep

  &DFT

      BASIS_SET_FILE_NAME BASIS_MOLOPT

      POTENTIAL_FILE_NAME GTH_POTENTIALS            

      CHARGE 0

      MULTIPLICITY 1

    &MGRID

       CUTOFF [Ry] 400 

    &END

    &QS

       METHOD GPW 

       EPS_DEFAULT 1.0E-6 

    &END

    &POISSON                      # POISSON solver for non-periodic calculation

       PERIODIC NONE

       PSOLVER WAVELET

    &END

    &SCF                              

      SCF_GUESS ATOMIC ! can be used to RESTART an interrupted calculation

      MAX_SCF 300

      EPS_SCF 1.0E-6 ! accuracy of the SCF procedure typically 1.0E-6 - 1.0E-7

      &OT

        PRECONDITIONER FULL_SINGLE_INVERSE

        MINIMIZER DIIS

      &END OT

      &MIXING

         METHOD BROYDEN_MIXING

         ALPHA  0.4

         BETA   1.5

         NBROYDEN 8

      &END MIXING

      &OUTER_SCF ! repeat the inner SCF cycle 10 times

        MAX_SCF 100

        EPS_SCF 1.0E-5 ! must match the above

        OPTIMIZER  DIIS

      &END

    &END SCF

    &XC

      &XC_FUNCTIONAL PBE

      &END XC_FUNCTIONAL 

    &END XC

  &END DFT

   &SUBSYS

    &CELL 

      ABC [angstrom] 10.00 10.00 10.00

    &END CELL

    &COORD

Na    0.0 0.0 0.0 

Cl    0.0 0.0 MYDIST

    &END COORD

    &TOPOLOGY

      CONNECTIVITY GENERATE

      &GENERATE

        BONDLENGTH_MAX 9

      &END

    &END

    &KIND Na                              

      BASIS_SET DZVP-MOLOPT-SR-GTH         

      POTENTIAL GTH-PBE          

    &END KIND

    &KIND Cl

      BASIS_SET DZVP-MOLOPT-SR-GTH

      POTENTIAL GTH-PBE

    &END KIND

  &END SUBSYS

&END FORCE_EVAL

The bash script used to get the calculation running is the following: 

rm -f ener_profile

#for dist in `seq 2.3 0.1 6.0`

for dist in 4.0 5.0 6.0 7.0 8.0 10.0; do

  #

  # compute energy 

  # 

  echo "Computing potential energy for distance $dist"

  sed "s/MYDIST/${dist}/g" mode1.inp > inp

  cp2k.popt inp > out 

  ener=`grep ' ENERGY| Total FORCE_EVAL' out | awk '{print $NF}'`

  echo $dist $ener >> ener_profile

done

[Patrick Gono]

Dear Mejdeddine,

You are using a non-periodic Poisson solver, suitable for this system of isolated atoms. However, the cell dimensions and atomic coordinates have to be chosen in such a way that the electronic density on the edges of the unit cell is negligible. Keep in mind that the cubic unit cell is situated with one corner in location (0, 0, 0) and the opposite corner in (10, 10, 10) Center the atoms in the middle of the cell, either manually, or by switching on the CENTER_COORDINATES statement in the TOPOLOGY section, which should center them in the middle of the unit cell by default.





Apart from that, I don't see any other immediate issue. What binding energy plot do you get? What is strange about it?





Yours sincerely,

Patrick Gono

--
You received this message because you are subscribed to the Google Groups "cp2k" group.
To unsubscribe from this group and stop receiving emails from it, send an email to cp2k+uns...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/cp2k/0798c015-5b94-4175-9179-71f3c2a5b5e3o%40googlegroups.com.

[mejdeddine mokhtar]

Thank's for your clarification. 

The curve should be showing a minimum of around 2.5 Angstrom and then tend to 1.5 eV even at a large distance.  What I found is completely wrong, even after setting the &TOPOLOGY &CENTER COOR.

I used the script below and the distance between both atoms is setting in the bash script.  Any help would be appreciated. 

The curve should look like the one I've attached below. 

To view this discussion on the web visit https://groups.google.com/d/msgid/cp2k/CAPL1_v95qaJijWmSD31i3%2B3RLU-QrDakeNe12fgb0DG6ELi%2BpQ%40mail.gmail.com.

[Lucas Lodeiro]

Hello Mejdeddine,

for non-periodic calculations, you need to set "PERIODIC NONE" in &CELL subsection too, to avoid pair interactions. Also, the cell need to be grosser than the system density. I guess at large distances the system will be Na⁺ and Cl⁻... for that the cell need to be, in the bond direction, the bond distance + twice the ionic radii of the biggest ion to ensure the density is zero at the borders, this is ~3.5 A due to the ionic radii of Cl⁻ is 1.7 A.

Using WAVELET solver needs a cubic box, so, if you want to compute the system at 10.0 A distance, the box will be 17.0A and cubic. With a cell of 10.0 A, you will have problems of non-zero density at 4.0 A distance, I guess.

Otherwise, the EPSs values must be equal for SCF and OUT_SCF, and the DEFAULT could be tightest (I use order of -6 and -12 respectively). And, the number of inner cycles, is useful less steps to use the outer cycles and improve the convergence.

Changing those things I get:

DIST (A) ENERGY(Ha)

1.0  -60.135198879741978
1.5  -61.855322129770521
2.0  -62.197077891951885
2.5  -62.227919856240661
3.0  -62.207241440696443
3.5  -62.182709877371700
4.0  -62.162755189608077
4.5  -62.148301846426023
5.0  -62.138459966222577
5.5  -62.131738730461414
6.0  -62.126565992226809
6.5  -62.122559283798338
7.0  -62.119778885873323

7.5  -62.118151081860589
8.0  -62.117267918233409

8.5  -62.116553100338109
9.0  -62.115268838231749

9.5  -62.114008578635961

10.0 -62.113128520139142

The minima is near 2.5 A, but I guess 2.6 A is closer... It is a coarse grid for the potential profile.

Also, the sentence "tend to 1.5 eV even at a large distance.", according to what I understand this value is arbitrary. At calculations this energy will be the energy of isolated Cl⁻ plus isolated Na⁺... But is a pseudopotential calculation, changing the pseudo or functional, the value will change. The important value is the energy depth of the minima with respect to the large converged distance. The large distance would be a coulomb r⁻¹ behavior.

I attach the input that I use to obtain list values.

Regards - Lucas Lodeiro

To view this discussion on the web visit https://groups.google.com/d/msgid/cp2k/CAHBAvwt1JXXfps0wB3CxPi7PpZxOcPn5BDbwjhrM204H8eoP5A%40mail.gmail.com.

[mejdeddine mokhtar]

I appreciate your answer. 

Unfornatullay, I used your values to plot the binding energy curve and I didn't get the right curve as it should. The bong length is around 2.5 however the minimum energy is wrong, It should be -4.2 eV.  

All the energies were referenced against the energy at a large distance (E - Eref) where Eref=E[-1].

The x-axis was the distances and the y-axis was (E -Eref) / Hartree. I've attached the plot. 

To view this discussion on the web visit https://groups.google.com/d/msgid/cp2k/CAOFT4PKcbV-0HiTsT7fo%3DU3WYbWm1YX4dUgiuqh33NLo3NZwgQ%40mail.gmail.com.

[mejdeddine mokhtar]

Hi Travis,

I meant that it should be -4.2 or -4.3, I found from another code (GPAW) the right curve but using cp2k I'm unable to get it. I think some parameters should be set correctly.  

Best.

Mejdeddine

Le mer. 8 juil. 2020 à 13:05, Travis <polla...@gmail.com> a écrit :

Hi,

I think you've converted incorrectly. Minimum at 2.5 Angstroms is -3.12 eV based on the data above.

-T

To view this discussion on the web visit https://groups.google.com/d/msgid/cp2k/eb1ccb3e-bbe7-4605-ac2b-eebdddafa6c0n%40googlegroups.com.

[Lucas Lodeiro]

The energy at 10.0 A is not a well behaved as "real large distance"... remember the ion interact throw coulomb interaction, charges of +1 and -1 have a interaction energy at 10.0 A of 1.4 eV, if we use this crude approximation, the energy of the minima is -4.5 eV respect to "converged large distance". Also you can calculate the isolated Na⁺ and Cl⁻, to get the infinite distance reference.

But if you want numeric accuracy, you need to use a better basis set at least, with diffuse and TZ, and maybe another functional.

To view this discussion on the web visit https://groups.google.com/d/msgid/cp2k/CAHBAvws6qYMiXEBgoHCWV9d3zbTt2b4%3Dt%3D0j4UYou-qgtp8g7Q%40mail.gmail.com.