Date Tags python

In order to do esp charge analysis, I generate charge density cube file and electric potential cube file from cp2k (with GPW method, GAPW method has a problem generate charge density cube file). I post here the steps I use horton and these cube files to do esp charge analysis (follow this link ).

Here is how I use it on killdevil.

module load anaconda
source activate my_root

Here is the input file for cp2k.

@SET RTYPE  MD
@SET RESTART  ON
@SET RESTARTFILE H2O64NaCl-1_500.restart
@SET MD_ENS NVT

&FORCE_EVAL
  METHOD QS
  &PROPERTIES
    #&RESP
    #  STRIDE 1 1 1
    #&END RESP
    &FIT_CHARGE
      FILENAME =NaCl_NVT_fit_charge
    &END FIT_CHARGE
  &END PROPERTIES
  &DFT
    &DENSITY_FITTING
    &END DENSITY_FITTING
    BASIS_SET_FILE_NAME ../../../cp2k-files/BASIS_MOLOPT
    POTENTIAL_FILE_NAME ../../../cp2k-files/GTH_POTENTIALS
    &MGRID
      NGRIDS 5
      CUTOFF 280
      REL_CUTOFF 40
    &END MGRID
    &SCF
      @IF ( ${RESTART} == OFF )
        SCF_GUESS ATOMIC
      @ENDIF
      @IF ( ${RESTART} == ON )
        SCF_GUESS RESTART
      @ENDIF
      &OT ON
        MINIMIZER DIIS
        ENERGY_GAP 0.001
        PRECONDITIONER FULL_ALL
      &END OT
      EPS_SCF      1.0E-6
      &PRINT
        &RESTART
        &END
      &END
    &END SCF
    &QS
      EXTRAPOLATION PS
      EXTRAPOLATION_ORDER 3
      METHOD GPW
    &END QS
    &XC
      &XC_FUNCTIONAL BLYP
      &END XC_FUNCTIONAL

      &XC_GRID
        XC_SMOOTH_RHO NN50
        XC_DERIV NN50_SMOOTH
      &END

      #&vdW_POTENTIAL
      #   POTENTIAL_TYPE PAIR_POTENTIAL
      #   &PAIR_POTENTIAL
      #      TYPE DFTD3
      #      PARAMETER_FILE_NAME ../../cp2k-files/dftd3.dat
      #      REFERENCE_FUNCTIONAL revPBE
      #   &END PAIR_POTENTIAL
      #&END vdW_POTENTIAL

    &END XC
    &PRINT
      &E_DENSITY_CUBE
        &EACH
          MD 100
        &END EACH
        STRIDE 1 1 1
      &END E_DENSITY_CUBE
      &V_HARTREE_CUBE
        &EACH
          MD 100
        &END EACH
        STRIDE 1 1 1
      &END V_HARTREE_CUBE
      &MULLIKEN
        &EACH
          MD 100
        &END EACH
      FILENAME =NaCl_NVT_mulliken
      &END MULLIKEN
      &LOWDIN
        &EACH
          MD 100
        &END EACH
      FILENAME =NaCl_NVT_lowdin
      &END LOWDIN
      #&HIRSCHFELD
      #  &EACH
      #    MD 100
      #  &END EACH
      #FILENAME =NaCl_NVT_hirschfeld
      #&END HIRSCHFELD
    &END PRINT
    #&LOCALIZE
    #  METHOD JACOBI
    #  #MAX_CRAZY_ANGLE 1.0
    #  MAX_ITER 100
    #  #EPS_LOCALIZATION 1E-8
    #  #MAX_CRAZY_ANGLE 0.05
    #  &PRINT
    #    &WANNIER_CENTERS
    #      &EACH
    #        MD 100
    #      &END EACH
    #      FILENAME =NaCl_NVT_wannier_centers.xyz
    #      IONS+CENTERS
    #    &END WANNIER_CENTERS
    #    &WANNIER_SPREADS
    #      &EACH
    #        MD 100
    #      &END EACH
    #      FILENAME =NaCl_NVT_wannier_spread.data
    #    &END WANNIER_SPREADS
    #  &END PRINT
    #  USE_HISTORY
    #&END LOCALIZE
  &END DFT
  &SUBSYS
    &CELL
      ABC 12.4934 12.4934 12.4934
    &END CELL
    &TOPOLOGY
      COORD_FILE_NAME  H2O64NaCl.xyz
      COORD_FILE_FORMAT xyz
    &END TOPOLOGY
    &COLVAR
     &DISTANCE
       ATOMS 193 194
     &END DISTANCE
     &PRINT
     &END
    &END COLVAR
    &KIND H
      BASIS_SET DZVP-MOLOPT-SR-GTH-q1
      POTENTIAL GTH-BLYP-q1
    &END KIND
    &KIND O
      BASIS_SET DZVP-MOLOPT-SR-GTH-q6
      POTENTIAL GTH-BLYP-q6
    &END KIND
    &KIND Na
      BASIS_SET DZVP-MOLOPT-SR-GTH-q9
      POTENTIAL GTH-BLYP-q9
    &END KIND
    &KIND Cl
      BASIS_SET DZVP-MOLOPT-SR-GTH-q7
      POTENTIAL GTH-BLYP-q7
    &END KIND
  &END SUBSYS
&END FORCE_EVAL

&GLOBAL
  PROJECT_NAME H2O64NaCl
  RUN_TYPE MD
  PRINT_LEVEL LOW

  &TIMINGS
     THRESHOLD 0.000001
  &END

  #WALLTIME 172500
  #WALLTIME 345000
  WALLTIME 1750

&END GLOBAL

&MOTION
  &CONSTRAINT
    &COLLECTIVE
      COLVAR 1
      INTERMOLECULAR
      TARGET [angstrom] 2.2
    &END COLLECTIVE
    &LAGRANGE_MULTIPLIERS
      COMMON_ITERATION_LEVELS 1
    &END
  &END CONSTRAINT
  &MD
    @IF ( ${MD_ENS} == NVE )
      ENSEMBLE NVE
    @ENDIF
    @IF ( ${MD_ENS} == NVT )
      ENSEMBLE NVT
      &THERMOSTAT
        REGION GLOBAL
        TYPE NOSE
        &NOSE
        &END NOSE
      &END THERMOSTAT
    @ENDIF
    STEPS 1
    TIMESTEP 0.5
    TEMPERATURE 298.0
  &END MD
  #&PRINT
  #    &TRAJECTORY  SILENT
  #    &EACH
  #      MD 1
  #    &END EACH
  #    FILENAME =NaCl_NVT.xyz
  #  &END TRAJECTORY
  #  &RESTART
  #    &EACH
  #      MD 1
  #    &END EACH
  #    ADD_LAST NUMERIC
  #  &END RESTART
  #  &VELOCITIES
  #    &EACH
  #      MD 1
  #    &END EACH
  #    FILENAME =NaCl_NVT_vel.xyz
  #  &END VELOCITIES
  #  &FORCES
  #    &EACH
  #      MD 1
  #    &END EACH
  #    FILENAME =NaCl_NVT_frc.xyz
  #  &END FORCES
  #&END PRINT
&END MOTION

&EXT_RESTART ON
  RESTART_DEFAULT F
  RESTART_FILE_NAME ${RESTARTFILE}
  RESTART_POS T
  RESTART_COUNTERS F
  @IF ( ${RTYPE} == MD )
    RESTART_VEL T
  @ENDIF
  @IF ( ${MD_ENS} == NVT )
    RESTART_THERMOSTAT T
  @ENDIF
&END EXT_RESTART

The cube files I have are H2O64NaCl-ELECTRON_DENSITY-1_0.cube and H2O64NaCl-v_hartree-1_0.cube.

This command is used to generate the cost function.

horton-esp-cost.py H2O64NaCl-v_hartree-1_0.cube cost.h5 --wdens=H2O64NaCl-ELECTRON_DENSITY-1_0.cube -s --wnear 8:2.0:2.0 1:2.0:2.0 11:2.0:2.0 17:2.0:2.0 --pbc=111 > esp-cost-out

It took me 6 hours to generate the cost function. After that, charges are fitted to the esp by this command.

horton-esp-fit.py cost.h5 charges.h5

Transfer it to csv file to read.

horton-hdf2csv.py charges.h5 charges.csv