Dear GTH,
I am preparing QM(DFT)/MM calculations for a chemical reaction catalyzed by an iron enzyme. I am interested in running the simulations both at BLYP and hybrid B3LYP level. While there is a Fe GTH optimized pseudopotential generated and available for the BLYP calculations in the CP2K database, there is no Fe basis set in the GTH_BASIS_SETS. Could you provide one? Can it be generated with the new ATOM BASIS_OPTIMIZATION codebase? Would you please address the same issue for B3LYP (BASIS/PSEUDOPOTENTIAL_OPTIMIZATION availability and accuracy)?
Thank you,
LC
Dear Juerg,
Thank you so much for your reply. I have one more question. Would you recommend using only MOLOPT basis sets for all elements (Fe, H, C, N, O) for the heme protein calculations rather than a mix of a MOLOPT basis set for Fe and more diffuse GTH basis sets for the other elements?
Best regards,
LaviniaDear Juerg,
I ran the regtest-admm for CH4 at B3LYP level (please see INPUT 1) successfully. Nonetheless, I get a get_gto_basis_set error (please see ERROR) when attempting to use ADMM for my system (please see INPUT 2). I am not sure what causes the error. Could you please provide insight? Thank you for all your assistance.
Sincerely,
Lavinia
ERROR:
***********************************************************
*** ERROR in get_gto_basis_set (MODULE basis_set_types) ***
***********************************************************
*** The pointer gto_basis_set is not associated ***
*** Program stopped at line number 433 of MODULE basis_set_types ***
===== Routine Calling Stack =====
4 hfx_create
3 quickstep_create_force_env
2 qmmm_create_force_env
1 CP2K
INPUT 1:
&FORCE_EVAL
METHOD Quickstep
&DFT
BASIS_SET_FILE_NAME ./BASIS_MOLOPT
POTENTIAL_FILE_NAME ./GTH_POTENTIALS
&MGRID
CUTOFF 100
REL_CUTOFF 30
&END MGRID
&QS
METHOD GPW
EPS_PGF_ORB 1.0E-12
EPS_FILTER_MATRIX 0.0e0
&END QS
&AUXILIARY_DENSITY_MATRIX_METHOD
METHOD BASIS_PROJECTION
ADMM_PURIFICATION_METHOD MO_DIAG
&END
&POISSON
PERIODIC NONE
PSOLVER MT
&END
&SCF
EPS_SCF 1.0E-6
SCF_GUESS ATOMIC
MAX_SCF 30
&OT ON
&END
&END SCF
&XC
&XC_FUNCTIONAL
&LYP
SCALE_C 0.81
&END
&BECKE88
SCALE_X 0.72
&END
&VWN
FUNCTIONAL_TYPE VWN3
SCALE_C 0.19
&END
&XALPHA
SCALE_X 0.08
&END
&END XC_FUNCTIONAL
&HF
&SCREENING
EPS_SCHWARZ 1.0E-10
SCREEN_ON_INITIAL_P FALSE
&END
&MEMORY
MAX_MEMORY 900
EPS_STORAGE_SCALING 0.1
&END
&INTERACTION_POTENTIAL
POTENTIAL_TYPE COULOMB
&END
FRACTION 0.20
&END
&END XC
&END DFT
&SUBSYS
&CELL
ABC 8.0 8.0 8.0
PERIODIC NONE
&END CELL
&COORD
C 0.0000 0.0000 0.0000
H 0.6297 0.6297 0.6297
H -0.6297 -0.6297 0.6297
H -0.6297 0.6297 -0.6297
H 0.6297 -0.6297 -0.6297
&END COORD
&KIND H
BASIS_SET DZVP-MOLOPT-SR-GTH-q1
AUX_FIT_BASIS_SET DZVP-MOLOPT-SR-GTH-q1
POTENTIAL GTH-BLYP-q1
&END KIND
&KIND C
BASIS_SET DZVP-MOLOPT-SR-GTH-q4
AUX_FIT_BASIS_SET DZVP-MOLOPT-SR-GTH-q4
POTENTIAL GTH-BLYP-q4
&END KIND
&END SUBSYS
&END FORCE_EVAL
&GLOBAL
PROJECT CH4-BP-MO_DIAG_B3LYP
PRINT_LEVEL LOW
RUN_TYPE MD
&TIMINGS
THRESHOLD 0.000000001
&END
&END GLOBAL
&MOTION
&MD
ENSEMBLE NVE
TIMESTEP 0.5
STEPS 2
&END
&END
INPUT 2:
METHOD GPW
# My scheme
EPS_DEFAULT 1.0E-12
EPS_PGF_ORB 1.0E-32
EPS_FILTER_MATRIX 0.0E+0
&END QS
&AUXILIARY_DENSITY_MATRIX_METHOD
METHOD BASIS_PROJECTION
ADMM_PURIFICATION_METHOD MO_DIAG
&END
SCREEN_ON_INITIAL_P FALSE
&END
&MEMORY
MAX_MEMORY 1300
EPS_STORAGE_SCALING 1.0E-1
&END
&INTERACTION_POTENTIAL
POTENTIAL_TYPE COULOMB
AUX_FIT_BASIS_SET DZVP-MOLOPT-SR-GTH-q16
&GLOBAL
PROJECT FeS_5
RUN_TYPE GEO_OPT
PRINT_LEVEL MEDIUM
WALLTIME 86400
&END GLOBAL
&FORCE_EVAL
METHOD Quickstep
&DFT
BASIS_SET_FILE_NAME /share/apps/cp2k/cp2k/tests/QS/BASIS_MOLOPT
POTENTIAL_FILE_NAME /share/apps/cp2k/cp2k/tests/QS/GTH_POTENTIALS
CHARGE 0
MULTIPLICITY 5
UKS T
&MGRID
CUTOFF 400
REL_CUTOFF 60
&END MGRID
&QS
METHOD GPW
EPS_DEFAULT 1.0E-12
&END QS
&SCF
MAX_SCF 500
&OUTER_SCF
EPS_SCF 1.0E-6
MAX_SCF 60
&END OUTER_SCF
&END SCF
&POISSON
PERIODIC NONE
PSOLVER WAVELET
&END POISSON
&XC
&XC_FUNCTIONAL
&LYP
SCALE_C 0.81
&END
&BECKE88
SCALE_X 0.72
&END
&VWN
FUNCTIONAL_TYPE VWN5
SCALE_C 0.19
&END
&END XC_FUNCTIONAL
&HF
&SCREENING
EPS_SCHWARZ 1.0E-10
&END
FRACTION 0.15
&END
&XC_GRID
XC_SMOOTH_RHO NN10
XC_DERIV SPLINE2_SMOOTH
&END XC_GRID
&END XC
&END DFT
&SUBSYS
&CELL
ABC 7.00 7.00 7.00
PERIODIC NONE
&END CELL
&COORD
Fe 0.0000 0.0000 0.0000
S 0.0000 0.0000 1.0600
&END COORD
&TOPOLOGY
&CENTER_COORDINATES
&END
&END
&KIND Fe
BASIS_SET DZVP-MOLOPT-SR-GTH
POTENTIAL GTH-PADE-q8
&END KIND
&KIND S
BASIS_SET DZVP-MOLOPT-GTH
POTENTIAL GTH-PADE-q6
&KIND Fe
BASIS_SET DZVP-MOLOPT-SR-GTH
POTENTIAL GTH-PADE-q8
&END KIND
&KIND S
BASIS_SET DZVP-MOLOPT-GTH
POTENTIAL GTH-PADE-q6
&END KIND
&KIND Fe
BASIS_SET DZVP-MOLOPT-SR-GTH
POTENTIAL GTH-BLYP-q16
&END KIND
&KIND S
BASIS_SET DZVP-MOLOPT-GTH
POTENTIAL GTH-BLYP-q6
&END KIND
gives me very different energy. From the PADE PP, the energy is: -28.4161411776336 while the one from BLYP gives the energy: -129.19258676213. Even the energy of the reaction Fe + S --> FeS from the BLYP PP is unphysical (4500 kJ/mol). I did not change anything else in the input file. Surprisingly this is not the case with the sulfur. For sulfur, with PADE PP the energy is: -9.889958274622 while with BLYP PP it is: -9.28358433022577. For Fe, with PADE PP it is: -18.4228633044091 and with BLYP PP it is: -121.648598758824
I wonder if it has something to do with the transition metals. Has anyone else experience similar problems. Will someone provide information on how to deal with transition metals in a QM calculation (non-periodic) in CP2K.
Thanks,
Megha