comp-chem/cp2k
1
0
Fork 0
mirror of https://github.com/cp2k/cp2k.git synced 2026-07-28 22:25:32 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

CDFT forces based on Hirshfeld partitioning of the electron density (#2111)

Browse source 
This update will allow for CDFT force calculations to be performed with Hirshfeld partitioning of the electron density, which provides a more robust description of atomic charges than Becke partitioning.
This commit is contained in:
chrisahart 2022-05-17 15:21:06 +01:00 • committed by GitHub
parent 6fcab5af1a
commit d19b00a1d0
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
14 changed files with 1218 additions and 445 deletions
Download patch file Download diff file Expand all files Collapse all files

43
src/input_cp2k_dft.F
View file

@ -4285,6 +4285,15 @@ CONTAINS
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="IN_MEMORY", &
description="Precompute gradients due to constraint during"// &
" initial formation of constraint and store them in memory. Does"// &
" nothing if forces are not calculated.", &
usage="IN_MEMORY", &
default_l_val=.FALSE., lone_keyword_l_val=.TRUE.)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
END SUBROUTINE create_cdft_control_section
! **************************************************************************************************
@ -4356,6 +4365,40 @@ CONTAINS
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="PRINT_DENSITY", &
description="Logical to control printing of Hirshfeld densities to .cube file.", &
usage="PRINT_DENSITY TRUE", &
default_l_val=.FALSE., lone_keyword_l_val=.TRUE.)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="ATOMS_MEMORY", &
description="Number of atomic gradients to store in memory.", &
usage="ATOMS_MEMORY", &
n_var=1, type_of_var=integer_t, &
default_i_val=80)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="USE_ATOMIC_CUTOFF", &
description="Logical to control use of ATOMIC_CUTOFF.", &
usage="USE_ATOMIC_CUTOFF TRUE", &
default_l_val=.TRUE., lone_keyword_l_val=.TRUE.)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="EPS_CUTOFF", &
description="Numerical cutoff for calculation of weight function.", &
usage="EPS_CUTOFF {real} ", default_r_val=1.0e-12_dp)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="ATOMIC_CUTOFF", &
description="Numerical cutoff for calculation of Hirshfeld densities.", &
usage="ATOMIC_CUTOFF {real} ", default_r_val=1.0e-12_dp)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
END SUBROUTINE create_hirshfeld_constraint_section
! **************************************************************************************************

1222
src/qs_cdft_methods.F
View file

File diff suppressed because it is too large Load diff

46
src/qs_cdft_types.F
View file

@ -18,6 +18,8 @@ MODULE qs_cdft_types
USE hirshfeld_types, ONLY: hirshfeld_type,&
release_hirshfeld_type
USE input_constants, ONLY: becke_cutoff_global,&
outer_scf_becke_constraint,&
outer_scf_hirshfeld_constraint,&
outer_scf_none,&
radius_single,&
shape_function_gaussian
@ -98,12 +100,17 @@ MODULE qs_cdft_types
!> \param radius Gaussian radius parameter
!> \param shape_function the constraint type: atomic density or single Gaussian
!> \param use_bohr determines whether to use angstrom or bohr units for the radii of Gaussians
!> \param use_atomic_cutoff Logical to control use of ATOMIC_CUTOFF
!> \param atomic_cutoff Numerical cutoff for calculation of Hirshfeld densities
!> \param atoms_memory Number of atomic gradients to store in memory
!> \param eps_cutoff Numerical cutoff for calculation of weight function
!> \param print_density Logical to control printing of Hirshfeld densities to .cube file
!> \param hirshfeld_env auxiliary type storing information about the Gaussians
! **************************************************************************************************
TYPE hirshfeld_constraint_type
INTEGER :: gaussian_shape, shape_function
LOGICAL :: use_bohr
REAL(KIND=dp) :: radius
INTEGER :: gaussian_shape, shape_function, atoms_memory
LOGICAL :: use_bohr, print_density, use_atomic_cutoff
REAL(KIND=dp) :: radius, eps_cutoff, atomic_cutoff
REAL(KIND=dp), DIMENSION(:), POINTER :: radii
TYPE(hirshfeld_type), POINTER :: hirshfeld_env
END TYPE hirshfeld_constraint_type
@ -154,6 +161,7 @@ MODULE qs_cdft_types
!> in a mixed_env that holds multiple CDFT states
!> \param first_iteration a flag to mark the first iteration for printing of additional data
!> \param print_weight logical which determines if CDFT weight functions should be saved to a file
!> \param in_memory logical which determines if the gradients of the Becke potential should be
!> \param is_constraint list of logicals which determines if an atom is included in a constraint group
!> \param strength Lagrangian multipliers of the constraints
!> \param target target values of the constraints
@ -194,8 +202,19 @@ MODULE qs_cdft_types
! Atomic gradients of the weight function at every grid point
REAL(KIND=dp), POINTER, &
DIMENSION(:, :, :, :) :: gradients
REAL(KIND=dp), POINTER, &
DIMENSION(:, :, :, :) :: gradients_x
REAL(KIND=dp), POINTER, &
DIMENSION(:, :, :, :) :: gradients_y
REAL(KIND=dp), POINTER, &
DIMENSION(:, :, :, :) :: gradients_z
! The weight function of this constraint group
TYPE(pw_p_type) :: weight
TYPE(pw_p_type), DIMENSION(:), POINTER :: hw_rho_atomic
TYPE(pw_p_type), DIMENSION(:), POINTER :: hw_rho_atomic_dr
TYPE(pw_p_type), DIMENSION(:), POINTER :: hw_rho_atomic_charge
TYPE(pw_p_type) :: hw_rho_total_constraint
TYPE(pw_p_type) :: hw_rho_total
END TYPE cdft_group_type
TYPE cdft_control_type
@ -214,7 +233,7 @@ MODULE qs_cdft_types
atomic_charges, fragment_density, &
fragments_integrated, flip_fragment(2), &
transfer_pot, external_control, &
first_iteration, print_weight
first_iteration, print_weight, in_memory
LOGICAL, POINTER, DIMENSION(:) :: is_constraint
REAL(KIND=dp), DIMENSION(:), POINTER :: strength, TARGET, value
REAL(KIND=dp), POINTER, &
@ -371,6 +390,7 @@ CONTAINS
cdft_control%should_purge = .FALSE.
cdft_control%purge_history = .FALSE.
cdft_control%calculate_metric = .FALSE.
cdft_control%in_memory = .FALSE.
cdft_control%purge_freq = 0
cdft_control%nbad_conv = 0
cdft_control%purge_offset = 0
@ -451,8 +471,17 @@ CONTAINS
DEALLOCATE (cdft_control%group(i)%coeff)
IF (ALLOCATED(cdft_control%group(i)%d_sum_const_dR)) &
DEALLOCATE (cdft_control%group(i)%d_sum_const_dR)
IF (ASSOCIATED(cdft_control%group(i)%gradients)) &
DEALLOCATE (cdft_control%group(i)%gradients)
IF (cdft_control%type == outer_scf_becke_constraint) THEN
IF (ASSOCIATED(cdft_control%group(i)%gradients)) &
DEALLOCATE (cdft_control%group(i)%gradients)
ELSE IF (cdft_control%type == outer_scf_hirshfeld_constraint) THEN
IF (ASSOCIATED(cdft_control%group(i)%gradients_x)) &
DEALLOCATE (cdft_control%group(i)%gradients_x)
IF (ASSOCIATED(cdft_control%group(i)%gradients_y)) &
DEALLOCATE (cdft_control%group(i)%gradients_y)
IF (ASSOCIATED(cdft_control%group(i)%gradients_z)) &
DEALLOCATE (cdft_control%group(i)%gradients_z)
END IF
IF (ASSOCIATED(cdft_control%group(i)%integrated)) &
DEALLOCATE (cdft_control%group(i)%integrated)
END DO
@ -518,8 +547,13 @@ CONTAINS
ALLOCATE (hirshfeld_control)
hirshfeld_control%use_bohr = .FALSE.
hirshfeld_control%print_density = .FALSE.
hirshfeld_control%use_atomic_cutoff = .TRUE.
hirshfeld_control%radius = 3.0_dp
hirshfeld_control%eps_cutoff = 1.0e-12_dp
hirshfeld_control%atomic_cutoff = 1.0e-12_dp
hirshfeld_control%shape_function = shape_function_gaussian
hirshfeld_control%atoms_memory = 80
hirshfeld_control%gaussian_shape = radius_single
NULLIFY (hirshfeld_control%hirshfeld_env)
NULLIFY (hirshfeld_control%radii)

115
src/qs_cdft_utils.F
View file

@ -86,7 +86,8 @@ MODULE qs_cdft_utils
! *** Public subroutines ***
PUBLIC :: becke_constraint_init, read_becke_section, read_cdft_control_section
PUBLIC :: hfun_scale, hirshfeld_constraint_init, cdft_constraint_print
PUBLIC :: hfun_scale, hirshfeld_constraint_init, cdft_constraint_print, &
cdft_print_hirshfeld_density, cdft_print_weight_function
CONTAINS
@ -564,7 +565,6 @@ CONTAINS
ALLOCATE (cdft_control%group(k)%atoms(natoms))
ALLOCATE (cdft_control%group(k)%coeff(natoms))
NULLIFY (cdft_control%group(k)%weight%pw)
NULLIFY (cdft_control%group(k)%gradients)
NULLIFY (cdft_control%group(k)%integrated)
tot_natoms = tot_natoms + natoms
! Now parse
@ -726,16 +726,19 @@ CONTAINS
TYPE(qs_control_type), INTENT(INOUT) :: qs_control
TYPE(section_vals_type), POINTER :: cdft_control_section
INTEGER :: k, nvar
LOGICAL :: exists
TYPE(cdft_control_type), POINTER :: cdft_control
TYPE(section_vals_type), POINTER :: becke_constraint_section, &
TYPE(section_vals_type), POINTER :: becke_constraint_section, group_section, &
hirshfeld_constraint_section, &
outer_scf_section, print_section
NULLIFY (outer_scf_section, hirshfeld_constraint_section, becke_constraint_section, &
print_section)
print_section, group_section)
cdft_control => qs_control%cdft_control
CPASSERT(ASSOCIATED(cdft_control))
group_section => section_vals_get_subs_vals(cdft_control_section, "ATOM_GROUP")
CALL section_vals_get(group_section, n_repetition=nvar, explicit=exists)
CALL section_vals_val_get(cdft_control_section, "TYPE_OF_CONSTRAINT", &
i_val=qs_control%cdft_control%type)
@ -771,11 +774,19 @@ CONTAINS
becke_constraint_section => section_vals_get_subs_vals(cdft_control_section, "BECKE_CONSTRAINT")
CALL section_vals_get(becke_constraint_section, explicit=exists)
IF (.NOT. exists) CPABORT("BECKE_CONSTRAINT section is missing.")
DO k = 1, nvar
NULLIFY (cdft_control%group(k)%gradients)
END DO
CALL read_becke_section(cdft_control, becke_constraint_section)
CASE (outer_scf_hirshfeld_constraint)
hirshfeld_constraint_section => section_vals_get_subs_vals(cdft_control_section, "HIRSHFELD_CONSTRAINT")
CALL section_vals_get(hirshfeld_constraint_section, explicit=exists)
IF (.NOT. exists) CPABORT("HIRSHFELD_CONSTRAINT section is missing.")
DO k = 1, nvar
NULLIFY (cdft_control%group(k)%gradients_x)
NULLIFY (cdft_control%group(k)%gradients_y)
NULLIFY (cdft_control%group(k)%gradients_z)
END DO
CALL read_hirshfeld_constraint_section(cdft_control, hirshfeld_constraint_section)
CASE DEFAULT
CPABORT("Unknown constraint type.")
@ -796,7 +807,6 @@ CONTAINS
!> \brief reads the input parameters needed for Hirshfeld constraint
!> \param cdft_control the cdft_control which holds the Hirshfeld constraint
!> \param hirshfeld_section the input section for a Hirshfeld constraint
!> \author Nico Holmberg [09.2018]
! **************************************************************************************************
SUBROUTINE read_hirshfeld_constraint_section(cdft_control, hirshfeld_section)
TYPE(cdft_control_type), INTENT(INOUT) :: cdft_control
@ -814,6 +824,11 @@ CONTAINS
CALL section_vals_val_get(hirshfeld_section, "GAUSSIAN_SHAPE", i_val=hirshfeld_control%gaussian_shape)
CALL section_vals_val_get(hirshfeld_section, "GAUSSIAN_RADIUS", r_val=hirshfeld_control%radius)
CALL section_vals_val_get(hirshfeld_section, "USE_BOHR", l_val=hirshfeld_control%use_bohr)
CALL section_vals_val_get(hirshfeld_section, "USE_ATOMIC_CUTOFF", l_val=hirshfeld_control%use_atomic_cutoff)
CALL section_vals_val_get(hirshfeld_section, "PRINT_DENSITY", l_val=hirshfeld_control%print_density)
CALL section_vals_val_get(hirshfeld_section, "EPS_CUTOFF", r_val=hirshfeld_control%eps_cutoff)
CALL section_vals_val_get(hirshfeld_section, "ATOMIC_CUTOFF", r_val=hirshfeld_control%atomic_cutoff)
IF (.NOT. hirshfeld_control%use_bohr) THEN
hirshfeld_control%radius = cp_unit_from_cp2k(hirshfeld_control%radius, "angstrom")
END IF
@ -837,19 +852,21 @@ CONTAINS
END SUBROUTINE read_hirshfeld_constraint_section
! **************************************************************************************************
! **************************************************************************************************
!> \brief Calculate fout = fun1/fun2 or fout = fun1*fun2
!> \param fout the output 3D potential
!> \param fun1 the first input 3D potential
!> \param fun2 the second input 3D potential
!> \param divide logical that decides whether to divide or multiply the input potentials
!> \param small customisable parameter to determine lower bound of division
! **************************************************************************************************
SUBROUTINE hfun_scale(fout, fun1, fun2, divide)
SUBROUTINE hfun_scale(fout, fun1, fun2, divide, small)
REAL(KIND=dp), DIMENSION(:, :, :), INTENT(OUT) :: fout
REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: fun1, fun2
LOGICAL, INTENT(IN) :: divide
REAL(KIND=dp), INTENT(IN) :: small
REAL(KIND=dp), PARAMETER :: small = 1.0e-12_dp
CHARACTER(len=*), PARAMETER :: routineN = 'hfun_scale', routineP = moduleN//':'//routineN
INTEGER :: i1, i2, i3, n1, n2, n3
@ -1250,9 +1267,7 @@ CONTAINS
! **************************************************************************************************
!> \brief Prints CDFT weight functions to cube files
!> \param qs_env the qs_env where to build the constraint
!> \par History
!> Created 9.2018 [Nico Holmberg]
!> \param qs_env ...
! **************************************************************************************************
SUBROUTINE cdft_print_weight_function(qs_env)
TYPE(qs_environment_type), POINTER :: qs_env
@ -1300,4 +1315,82 @@ CONTAINS
END SUBROUTINE cdft_print_weight_function
! **************************************************************************************************
!> \brief Prints Hirshfeld weight function and promolecule density
!> \param qs_env ...
! **************************************************************************************************
SUBROUTINE cdft_print_hirshfeld_density(qs_env)
TYPE(qs_environment_type), POINTER :: qs_env
CHARACTER(len=*), PARAMETER :: routineN = 'cdft_print_hirshfeld_density', &
routineP = moduleN//':'//routineN
CHARACTER(LEN=default_path_length) :: middle_name
INTEGER :: iatom, igroup, unit_nr
LOGICAL :: mpi_io
TYPE(cdft_control_type), POINTER :: cdft_control
TYPE(cp_logger_type), POINTER :: logger
TYPE(cp_para_env_type), POINTER :: para_env
TYPE(dft_control_type), POINTER :: dft_control
TYPE(particle_list_type), POINTER :: particles
TYPE(pw_env_type), POINTER :: pw_env
TYPE(qs_subsys_type), POINTER :: subsys
TYPE(section_vals_type), POINTER :: cdft_constraint_section
NULLIFY (cdft_constraint_section, logger, particles, dft_control, &
para_env, subsys, cdft_control, pw_env)
logger => cp_get_default_logger()
CALL get_qs_env(qs_env, subsys=subsys, para_env=para_env, dft_control=dft_control, pw_env=pw_env)
CALL qs_subsys_get(subsys, particles=particles)
cdft_control => dft_control%qs_control%cdft_control
cdft_constraint_section => section_vals_get_subs_vals(qs_env%input, "DFT%QS%CDFT")
mpi_io = .TRUE.
DO igroup = 1, SIZE(cdft_control%group)
middle_name = "hw_rho_total"//TRIM(ADJUSTL(cp_to_string(igroup)))
unit_nr = cp_print_key_unit_nr(logger, cdft_constraint_section, "PROGRAM_RUN_INFO", mpi_io=mpi_io, &
file_position="REWIND", middle_name=middle_name, extension=".cube")
CALL cp_pw_to_cube(cdft_control%group(igroup)%hw_rho_total%pw, unit_nr, "CDFT Weight Function", mpi_io=mpi_io, &
particles=particles, stride=section_get_ivals(cdft_constraint_section, "PROGRAM_RUN_INFO%WEIGHT_FUNCTION%STRIDE"))
CALL cp_print_key_finished_output(unit_nr, logger, cdft_constraint_section, "PROGRAM_RUN_INFO", mpi_io=mpi_io)
END DO
DO igroup = 1, SIZE(cdft_control%group)
middle_name = "hw_rho_total_constraint_"//TRIM(ADJUSTL(cp_to_string(igroup)))
unit_nr = cp_print_key_unit_nr(logger, cdft_constraint_section, "PROGRAM_RUN_INFO", mpi_io=mpi_io, &
file_position="REWIND", middle_name=middle_name, extension=".cube")
CALL cp_pw_to_cube(cdft_control%group(igroup)%hw_rho_total_constraint%pw, unit_nr, &
"CDFT Weight Function", mpi_io=mpi_io, particles=particles, &
stride=section_get_ivals(cdft_constraint_section, "PROGRAM_RUN_INFO%WEIGHT_FUNCTION%STRIDE"))
CALL cp_print_key_finished_output(unit_nr, logger, cdft_constraint_section, "PROGRAM_RUN_INFO", mpi_io=mpi_io)
END DO
DO igroup = 1, SIZE(cdft_control%group)
DO iatom = 1, (cdft_control%natoms)
middle_name = "hw_rho_atomic_"//TRIM(ADJUSTL(cp_to_string(iatom)))
unit_nr = cp_print_key_unit_nr(logger, cdft_constraint_section, "PROGRAM_RUN_INFO", mpi_io=mpi_io, &
file_position="REWIND", middle_name=middle_name, extension=".cube")
CALL cp_pw_to_cube(cdft_control%group(igroup)%hw_rho_atomic(iatom)%pw, unit_nr, &
"CDFT Weight Function", mpi_io=mpi_io, particles=particles, &
stride=section_get_ivals(cdft_constraint_section, "PROGRAM_RUN_INFO%WEIGHT_FUNCTION%STRIDE"))
CALL cp_print_key_finished_output(unit_nr, logger, cdft_constraint_section, "PROGRAM_RUN_INFO", mpi_io=mpi_io)
END DO
END DO
END SUBROUTINE cdft_print_hirshfeld_density
END MODULE qs_cdft_utils

13
tests/QS/regtest-cdft-hirshfeld-2/TEST_FILES
View file

@ -3,16 +3,15 @@
# see regtest/TEST_FILES
#
HeH-noconstraint.inp 1 2e-13 -3.01067446615063
HeH-cdft-state-1.inp 71 5e-08 0.161058337228
HeH-cdft-state-2.inp 71 1e-09 1.849710970513
HeH-cdft-state-1.inp 71 1e-7 0.161058337228
HeH-cdft-state-2.inp 71 1e-7 1.849710970513
# These tests compute the electronic coupling and related quantities
# Coupling via rotation
HeH-mixed-cdft-1.inp 73 5e-10 204.357138287654
HeH-mixed-cdft-1.inp 73 1e-7 204.357138287654
# Coupling via Lowdin orthogonalization
HeH-mixed-cdft-2.inp 76 5e-10 154.149819062673
HeH-mixed-cdft-2.inp 76 1e-7 154.149819062673
# CDFT-CI
HeH-mixed-cdft-3.inp 77 3e-12 -2.80477609770094
HeH-mixed-cdft-3.inp 77 1e-7 -2.80477609770094
# Rerun some tests using parallel mode (NGROUPS 2)
# Electronic coupling and related
HeH-mixed-cdft-8.inp 73 5e-10 204.357138287651
#EOF
HeH-mixed-cdft-8.inp 73 1e-7 204.357138287651

32
tests/QS/regtest-cdft-hirshfeld-3/HeH-cdft-1.inp Normal file
View file

@ -0,0 +1,32 @@
@SET RESTART_WFN TRUE
@SET WFN_FILE HeH-noconstraint-1_0.wfn
@SET PROJECT_NAME HeH-cdft-1
@SET WRITE_WFN 0
@SET CHARGE 1
@SET WRITE_CUBE FALSE
@SET CENTER_SYS ON
@SET XYZFILE HeH.xyz
@SET HIRSHFELD_ACTIVE TRUE
@SET HIRSHFELD_FRAGMENT FALSE
@SET PRINT_WEIGHT_FUNCTION FALSE
@SET MAX_SCF 0
@SET HIRSHFELD_TARGET 1.0
@SET HIRSHFELD_STR 0.186894372937
@SET HIRSHFELD_SHAPEFN DENSITY
&GLOBAL
PROJECT ${PROJECT_NAME}
RUN_TYPE ENERGY_FORCE
PRINT_LEVEL MEDIUM
PREFERRED_DIAG_LIBRARY SL
&END GLOBAL
&FORCE_EVAL
METHOD QS
@SET HIRSHFELD_STR ${HIRSHFELD_STR}
@SET HIRSHFELD_TARGET ${HIRSHFELD_TARGET}
@include dft-common-params.inc
@include subsys.inc
&END FORCE_EVAL

23
tests/QS/regtest-cdft-hirshfeld-3/HeH-noconstraint.inp Normal file
View file

@ -0,0 +1,23 @@
@SET RESTART_WFN FALSE
@SET WFN_FILE HeH-noconstraint-1_0.wfn
@SET PROJECT_NAME HeH-noconstraint
@SET WRITE_WFN 1
@SET CHARGE 1
@SET WRITE_CUBE FALSE
@SET CENTER_SYS ON
@SET XYZFILE HeH.xyz
@SET HIRSHFELD_ACTIVE FALSE
&GLOBAL
PROJECT ${PROJECT_NAME}
RUN_TYPE ENERGY_FORCE
PRINT_LEVEL MEDIUM
PREFERRED_DIAG_LIBRARY SL
&END GLOBAL
&FORCE_EVAL
METHOD QS
@include dft-common-params.inc
@include subsys.inc
&END FORCE_EVAL

4
tests/QS/regtest-cdft-hirshfeld-3/HeH.xyz Normal file
View file

@ -0,0 +1,4 @@
2
He 3.000000 3.000000 2.600000
H 3.000000 3.000000 3.400000

7
tests/QS/regtest-cdft-hirshfeld-3/TEST_FILES Normal file
View file

@ -0,0 +1,7 @@
# runs are executed in the same order as in this file
# the second field tells which test should be run in order to compare with the last available output
# see regtest/TEST_FILES
#
HeH-noconstraint.inp 72 1e-7 0.1450972684448
HeH-cdft-1.inp 72 1e-7 0.1552195046628
#EOF

82
tests/QS/regtest-cdft-hirshfeld-3/dft-common-params.inc Normal file
View file

@ -0,0 +1,82 @@
&DFT
@IF ( ${HIRSHFELD_ACTIVE} == TRUE )
@include hirshfeld_qs.inc
@ENDIF
@IF ( ${HIRSHFELD_ACTIVE} == FALSE )
&QS
METHOD GPW
EPS_DEFAULT 1.0E-12
EXTRAPOLATION ASPC
EXTRAPOLATION_ORDER 3
&END QS
@ENDIF
BASIS_SET_FILE_NAME BASIS_MOLOPT
POTENTIAL_FILE_NAME POTENTIAL
@IF ( ${RESTART_WFN} == TRUE )
WFN_RESTART_FILE_NAME ${WFN_FILE}
@ENDIF
LSD
CHARGE ${CHARGE}
&MGRID
CUTOFF 100
NGRIDS 5
&END MGRID
&SCF
@IF ( ${RESTART_WFN} == TRUE )
SCF_GUESS RESTART
@ENDIF
@IF ( ${RESTART_WFN} == FALSE )
SCF_GUESS ATOMIC
@ENDIF
EPS_SCF 1.0E-5
CHOLESKY INVERSE_DBCSR
MAX_SCF 20
&OT ON
! DIIS convergence might not be constant with different mpiranks => use CG
MINIMIZER CG
PRECONDITIONER FULL_ALL
ALGORITHM IRAC
&END OT
&OUTER_SCF ON
EPS_SCF 1.0E-5
MAX_SCF 2
&END
&PRINT
&RESTART
FILENAME ./${PROJECT_NAME}
BACKUP_COPIES 0
COMMON_ITERATION_LEVELS 1
&EACH
JUST_ENERGY ${WRITE_WFN}
QS_SCF 0
&END EACH
&END RESTART
&RESTART_HISTORY OFF
&END RESTART_HISTORY
&END PRINT
&END SCF
&XC
&XC_FUNCTIONAL PBE
&END XC_FUNCTIONAL
&XC_GRID
XC_DERIV SPLINE2
XC_SMOOTH_RHO NONE
&END XC_GRID
&END XC
&PRINT
&MULLIKEN OFF
&END MULLIKEN
&HIRSHFELD OFF
&END HIRSHFELD
@IF ( ${WRITE_CUBE} == TRUE )
&E_DENSITY_CUBE ON
STRIDE 1 1 1
&END E_DENSITY_CUBE
@ENDIF
&END PRINT
&END DFT
&PRINT
&FORCES ON
NDIGITS 13
&END
&END

46
tests/QS/regtest-cdft-hirshfeld-3/hirshfeld_qs.inc Normal file
View file

@ -0,0 +1,46 @@
&QS
METHOD GPW
EPS_DEFAULT 1.0E-12
EXTRAPOLATION ASPC
EXTRAPOLATION_ORDER 3
&CDFT
TYPE_OF_CONSTRAINT HIRSHFELD
STRENGTH ${HIRSHFELD_STR}
TARGET ${HIRSHFELD_TARGET}
&ATOM_GROUP
ATOMS 1 2
COEFF 1 -1
&END ATOM_GROUP
&OUTER_SCF ON
TYPE CDFT_CONSTRAINT
MAX_SCF 50
EPS_SCF 1.0E-3
OPTIMIZER NEWTON_LS
&CDFT_OPT ON
CONTINUE_LS TRUE
&END CDFT_OPT
&END
&HIRSHFELD_CONSTRAINT
SHAPE_FUNCTION DENSITY
&END HIRSHFELD_CONSTRAINT
&PROGRAM_RUN_INFO ON
! Print weight function to cube file
@IF ( ${PRINT_WEIGHT_FUNCTION} == TRUE )
&WEIGHT_FUNCTION
&END WEIGHT_FUNCTION
@ENDIF
&EACH
QS_SCF 1
&END EACH
COMMON_ITERATION_LEVELS 2
ADD_LAST NUMERIC
FILENAME ./${PROJECT_NAME}
&END PROGRAM_RUN_INFO
&END CDFT
&END QS

20
tests/QS/regtest-cdft-hirshfeld-3/subsys.inc Normal file
View file

@ -0,0 +1,20 @@
&SUBSYS
&CELL
PERIODIC XYZ
ABC 6. 6. 6.
&END CELL
&TOPOLOGY
COORD_FILE_FORMAT XYZ
COORD_FILE_NAME ${XYZFILE}
&CENTER_COORDINATES ${CENTER_SYS}
&END CENTER_COORDINATES
&END TOPOLOGY
&KIND He
BASIS_SET SZV-MOLOPT-SR-GTH
POTENTIAL GTH-PBE-q2
&END KIND
&KIND H
BASIS_SET SZV-MOLOPT-SR-GTH
POTENTIAL GTH-PBE-q1
&END KIND
&END SUBSYS

9
tests/QS/regtest-cdft-hirshfeld/TEST_FILES
View file

@ -6,8 +6,7 @@ HeH-noconstraint.inp 1 2e-13
He+-noconstraint.inp 1 2e-13 -2.05007934458372
H-noconstraint.inp 1 4e-13 -0.45734465780293
# These tests use different constraint formalisms so their value differs (see outputted charges)
HeH-cdft-1.inp 71 3e-11 1.135281643193
HeH-cdft-2.inp 71 2e-11 1.681172869379
HeH-cdft-3.inp 71 8e-08 1.197176189824
HeH-cdft-4.inp 71 2e-11 1.567706385507
#EOF
HeH-cdft-1.inp 71 1e-7 1.135281643193
HeH-cdft-2.inp 71 1e-7 1.681172869379
HeH-cdft-3.inp 71 1e-7 1.197176189824
HeH-cdft-4.inp 71 1e-7 1.567706385507

1
tests/TEST_DIRS
View file

@ -44,6 +44,7 @@ QS/regtest-negf
QS/regtest-negf-fft fftw3
QS/regtest-ri-rpa-axk libint
QS/regtest-ri-rpa-rse libint
QS/regtest-cdft-hirshfeld-3
QS/regtest-cdft-hirshfeld-2 parallel mpiranks>1
QS/regtest-cdft-hirshfeld
SIRIUS/regtest-1 sirius