This is an old revision of the document!

Validation of a KCl QMMM model

(exercise by Matthew Watkins, University college, London)

In this exercise you will validate the mixed quamtum/classical model for a KCl slab. The present exercise is referring to the following paper: 10.1002/jcc.23904.

You will optimize the geometry of a KCL slab with the same arrangement as the one depicted below

For simplicity, you will consider three layers.
You will compare the calculation performed with the full QM, one layer QM and two MM, two layer QM and one MM.
In particular you will compare the band gaps and the density of states.

1. Task: Prepare the input files

The file input.inp contains the partition of the QM and MM regions by specifying the atom index. Check in the kcl.xyz how the indexes are distributed. The slab orthogonal coordinate is the y.

Make three copies of input.inp and call them qm_1l.inp, qm_2l.inp and qm.inp.
qm_1l.inp should have PROJECT KCl_1 and one QM layer. Modify accordingly by looking first at the kcl.xyz coordinate files.
qm_2l.inp should have PROJECT KCl_2 and two QM layers. Modify accordingly by looking first at the kcl.xyz coordinate files.
qm.inp should have PROJECT KCl_QM and the full QM treatment. For this, it is sufficient to change QMMM at the beginning to QS.

Questions

Required Files

When you are dealing with complex job structure, the input can be simplified by splitting it into multiple files. We are going to use separate files for the coordinates, the QM part, the MM part. All these files should reside in the same directory as the main input file.

The provided files are all in the directory /home/psd/Exercise_13. When you create new input files with different parameters, remember to change the name of the PROJECT as well.

input.inp

@SET METHOD = QMMM # FIST all classical treatment # QS all quantum treatment
 
&GLOBAL
  FLUSH_SHOULD_FLUSH
  PRINT_LEVEL LOW
  PROJECT KCl
  RUN_TYPE GEO_OPT
&END GLOBAL
 
&FORCE_EVAL
  METHOD $METHOD
  @include QS.inc
  @include MM.inc
  &QMMM
    #this defines the QS cell in the QMMM calc
    &CELL
      ABC 12.6 15.0 12.6
      PERIODIC XZ
    &END CELL
    ECOUPL GAUSS # use GEEP method
    NOCOMPATIBILITY
    USE_GEEP_LIB 6  # use GEEP method
    &PERIODIC # apply periodic potential
      #in this case QM box = MM box in XZ so turn
      #off coupling/recoupling of the QM multipole
      &MULTIPOLE OFF
      &END
    &END PERIODIC
    #these are just the ionic radii of K Cl
    #but should be treated as parameters in general
    #fit to some physical property
    &MM_KIND K
      RADIUS 1.52
    &END MM_KIND
    &MM_KIND Cl
      RADIUS 1.67
    &END MM_KIND
    #define the model
    &QM_KIND K
      MM_INDEX 25..32 41..48
    &END QM_KIND
        &MM_KIND Cl
      RADIUS 1.67
    &END MM_KIND
    #define the model
    &QM_KIND K
      MM_INDEX 25..32 41..48
    &END QM_KIND
    &QM_KIND Cl
      MM_INDEX 17..24 33..40
    &END QM_KIND
  &END QMMM
 
  &SUBSYS
    #this defines the cell of the whole system
    #must be orthorhombic, I think
    &CELL
      ABC 12.6 100.0 12.6
    &END CELL
    &TOPOLOGY
      COORD_FILE_NAME kcl.xyz
      COORD_FILE_FORMAT XYZ
      &GENERATE
         &ISOLATED_ATOMS
         #ignores bonds dihedrals etc in classical part
            LIST 1..48
         &END
      &END
    &END
    &KIND K
      ELEMENT K
      BASIS_SET DZVP-MOLOPT-SR-GTH
      POTENTIAL GTH-PBE-q9
    &END KIND
    &KIND Cl
      BASIS_SET DZVP-MOLOPT-GTH
      POTENTIAL GTH-PBE-q7
    &END
  &END SUBSYS
&END FORCE_EVAL
 
#should be able to use most motion sections
#analytic stress tensor not available, I think
@include motion.inc

and includes as separate files, using the @include macro, for the QS, MM and motion sections

QS.inc

  &DFT
    BASIS_SET_FILE_NAME BASIS_MOLOPT
    POTENTIAL_FILE_NAME GTH_POTENTIALS
    &MGRID
      COMMENSURATE
      CUTOFF 150
    &END MGRID
    &QS
      EPS_DEFAULT 1.0E-12
    &END QS
    &SCF
      EPS_SCF 1.0E-06
      MAX_SCF 26
      SCF_GUESS RESTART
      &OT
        MINIMIZER CG
        PRECONDITIONER FULL_SINGLE_INVERSE
        ENERGY_GAP 0.001
      &END OT
      &OUTER_SCF
        EPS_SCF 1.0E-05
      &END OUTER_SCF
    &END SCF
    &XC
      &XC_FUNCTIONAL PBE
      &END XC_FUNCTIONAL
    &END XC
    &PRINT
       &MO_CUBES
           NLUMO 10
           WRITE_CUBE T
       &END MO_CUBES
       &V_HARTREE_CUBE
           STRIDE 2 2 2
       &END
    &END PRINT
  &END DFT

MM.inc

  &MM
    &FORCEFIELD
      &CHARGE
         ATOM K
         CHARGE 1.0
      &END CHARGE
      &CHARGE
         ATOM Cl
         CHARGE -1.0
      &END CHARGE
      &NONBONDED
        &WILLIAMS
          atoms K   Cl
          A [eV] 4117.9
          B [angstrom^-1] 3.2808
          C [eV*angstrom^6] 0.0
          RCUT [angstrom] 3.0
        &END WILLIAMS
        &WILLIAMS
          atoms Cl  Cl
          A [eV] 1227.2
          B [angstrom^-1] 3.1114
          C [eV*angstrom^6] 124.0
          RCUT [angstrom] 3.0
        &END WILLIAMS
        &WILLIAMS
          atoms K   K
          A [eV] 3796.9
          B [angstrom^-1] 3.84172
          C [eV*angstrom^6] 124.0
          RCUT [angstrom] 3.0
        &END WILLIAMS
      &END NONBONDED
    &END FORCEFIELD
    &POISSON
      &EWALD
        EWALD_TYPE spme
        ALPHA .44
        GMAX  40
      &END EWALD
    &END POISSON
  &END MM

motion.inc

&MOTION
  &GEO_OPT
     OPTIMIZER LBFGS
  &END
  &CONSTRAINT
     &FIXED_ATOMS
     LIST 1..16
     EXCLUDE_MM .FALSE.
     EXCLUDE_QM .TRUE.
     &END FIXED_ATOMS
  &END CONSTRAINT
&END MOTION

kcl.xyz

48
 
Cl         0.00000       15.00000        0.00000
Cl         3.15000       15.00000        3.15000
Cl         0.00000       15.00000        6.30000
Cl         3.15000       15.00000        9.45000
Cl         6.30000       15.00000        0.00000
Cl         9.45000       15.00000        3.15000
Cl         6.30000       15.00000        6.30000
Cl         9.45000       15.00000        9.45000
K          3.15000       15.00000        0.00000
K          0.00000       15.00000        3.15000
K          3.15000       15.00000        6.30000
K          0.00000       15.00000        9.45000
K          9.45000       15.00000        0.00000
K          6.30000       15.00000        3.15000
K          9.45000       15.00000        6.30000
K          6.30000       15.00000        9.45000
Cl         3.15000       18.15000        0.00000
Cl         0.00000       18.15000        3.15000
Cl         3.15000       18.15000        6.30000
Cl         0.00000       18.15000        9.45000
Cl         9.45000       18.15000        0.00000
Cl         6.30000       18.15000        3.15000
Cl         9.45000       18.15000        6.30000
Cl         6.30000       18.15000        9.45000
K          0.00000       18.15000        0.00000
K          3.15000       18.15000        3.15000
K          0.00000       18.15000        6.30000
K          3.15000       18.15000        9.45000
K          6.30000       18.15000        0.00000
K          9.45000       18.15000        3.15000
K          6.30000       18.15000        6.30000
K          9.45000       18.15000        9.45000
Cl         0.00000       21.30000        0.00000
Cl         3.15000       21.30000        3.15000
Cl         0.00000       21.30000        6.30000
Cl         3.15000       21.30000        9.45000
Cl         6.30000       21.30000        0.00000
Cl         9.45000       21.30000        3.15000
Cl         6.30000       21.30000        6.30000
Cl         9.45000       21.30000        9.45000
K          3.15000       21.30000        0.00000
K          0.00000       21.30000        3.15000
K          3.15000       21.30000        6.30000
K          0.00000       21.30000        9.45000
K          9.45000       21.30000        0.00000
K          6.30000       21.30000        3.15000
K          9.45000       21.30000        6.30000
K          6.30000       21.30000        9.45000