[CP2K-user] [CP2K:13546] IC-QMMM with single charge in front of a metal plane: open boundary corrections

[Katharina Doblhoff-Dier]

Dear Dorothea,
Thank you for your answer. 

I am not sure if I understand your question correctly 
>
 
On the longer run, we are interested in describing an electrolyte above a 
metal surface and using a charge imbalance to create a potential drop at 
the interface very much as one would do in a pure DFT calculation, but then 
withouth the need to describe the metal quantum mechanically. I know that 
this is not the original idea behind the image charge method as implemented 
in cp2k, but I thought it may offer a work around. So to understand how 
things work, I simply considered an H2+ molecule above a surface. This is 
similar to the model system used by Siepmann and Sprik in their 1995 
publication, just that they used an MM point charge above a metal. In 
principle, it should be possible to determine V_0 from the condition that 
the energy should be minimized, and the resulting counter charge should 
exactly cancel the charge in front of the metal.

and what your computational setup is.
>
I paste here an example of the input that I used with MT boundary 
conditions. The resulting plot for the Hartree-potential (obtained using  
&V_HARTREE_CUBE and then averaging over x and y to get a plot over z) is 
shown in the third image in my original post (green line). 
As you can see, I simply adapted the input file from the GUA example on the 
web, by replacing the GUA in the xyz file by H2, setting the DFT charge to 
1 and changing the poisson solver in the DFT and MM part. The potential 
EXT_POTENTIAL was determined by running 2 calculations with changing 
EXT_POTENTIAL_1 and EXT_POTENTIAL_2, and then determining EXT_POTENTIAL_3 
by requiring that the resulting charge should be -1 and assuming a linear 
relationship (which for this system is perfectly followed, giving me a 
total charge on the metal of -1)

-----------------Begin sample input file using the MT poisson 
solver------------------
&FORCE_EVAL
  METHOD QMMM
  &DFT
    CHARGE 1
    LSD
    BASIS_SET_FILE_NAME BASIS_MOLOPT
    POTENTIAL_FILE_NAME GTH_POTENTIALS
    &POISSON
      &EWALD
        EWALD_TYPE ewald
        ALPHA .44
        GMAX 21
      &END EWALD
      POISSON_SOLVER MT
      PERIODIC XY
    &END POISSON
    &MGRID
      COMMENSURATE
      CUTOFF 300
      NGRIDS 5
    &END MGRID
    &QS
      METHOD GPW
      EXTRAPOLATION ASPC
      EXTRAPOLATION_ORDER 3
    &END QS
    &SCF
      MAX_SCF 300
       SCF_GUESS ATOMIC
     # SCF_GUESS RESTART
      EPS_SCF 1.0E-6
      &OT
        PRECONDITIONER  FULL_SINGLE_INVERSE
        MINIMIZER DIIS
      &END
    &END SCF
    &XC
      &XC_FUNCTIONAL PBE
      &END XC_FUNCTIONAL
        &VDW_POTENTIAL
          DISPERSION_FUNCTIONAL PAIR_POTENTIAL
          &PAIR_POTENTIAL
           TYPE DFTD3
           CALCULATE_C9_TERM .TRUE.
           REFERENCE_C9_TERM
           PARAMETER_FILE_NAME dftd3.dat
           REFERENCE_FUNCTIONAL PBE
           R_CUTOFF [angstrom] 16.0
          &END PAIR_POTENTIAL
        &END VDW_POTENTIAL
    &END XC
    &PRINT
     &V_HARTREE_CUBE
     &END
    &END
  &END DFT

  &MM
    &FORCEFIELD
      &CHARGE
        ATOM Au
        CHARGE 0
      &END CHARGE
      &CHARGE
        ATOM H
        CHARGE 0
      &END CHARGE

      &SPLINE
       EPS_SPLINE 1.E-5
       #EMAX_SPLINE 2.0
      &END

       &NONBONDED
        &EAM
          atoms Au Au
          PARM_FILE_NAME Au.pot
        &END EAM


       &LENNARD-JONES
          atoms Au H
          EPSILON 0.0
          SIGMA 3.166
          RCUT  15
       &END LENNARD-JONES

       &LENNARD-JONES
         atoms H H
         EPSILON 0.0
         SIGMA 3.166
         RCUT 15
       &END LENNARD-JONES

      &END

    &END FORCEFIELD

    &POISSON
      &EWALD
        EWALD_TYPE ewald
        ALPHA .44
        GMAX 21
      &END EWALD
      POISSON_SOLVER MT
      PERIODIC XY
    &END POISSON
  &END MM

  &QMMM
    CENTER NEVER
    &CELL
      ABC [angstrom]   34.6055 29.9693 60.0
      PERIODIC XY
    &END CELL

    &QM_KIND H
     MM_INDEX 577..578
    &END QM_KIND

    &IMAGE_CHARGE
     EXT_POTENTIAL 1.24534
     MM_ATOM_LIST 1..576
     WIDTH 3.5
    &END IMAGE_CHARGE

    &PRINT
     &IMAGE_CHARGE_INFO
     &END
    &END

  &END QMMM

  &SUBSYS
    &CELL
      ABC [angstrom]   34.6055 29.9693 60.0
      PERIODIC XY
    &END CELL
    &TOPOLOGY
      COORD_FILE_NAME Au_gua_image_dampFunc-opt.xyz
      COORDINATE xyz
    &END
    &KIND H
      BASIS_SET DZVP-MOLOPT-SR-GTH
      POTENTIAL GTH-PBE-q1
    &END KIND

  &END SUBSYS
&END FORCE_EVAL

&GLOBAL
  PROJECT Au_gua_image_dampFunc
  RUN_TYPE energy
&END GLOBAL

-----------------End sample input file using the MT poisson 
solver------------------

>If you have a charged system in the QM part, then you need a solver like 
MT, like in a standard QM calculation.
Yes, I would have thought, that this is what I am doing, but the resulting 
plot for the hartree potential looks strange (shielding charges in the 
metal, leading to a potential gradient in the metal instead of a constant) 
and I cannot seem to find an energy minimum, when the potential is chosen 
such, that the image charge is exatly shielded, which was one of the test 
systems in the Siepmann Sprik paper 
(https://aip.scitation.org/doi/abstract/10.1063/1.469429).

I guess the latter is my main issue, as we will depend heavily on the 
energies. (Hence I need to understand what they contain)

Thank you for your help,
Katharina



-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.cp2k.org/archives/cp2k-user/attachments/20200624/519c671e/attachment.htm>