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