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

serial mode treatment of mixed CDFT states, more optimizer types, and

Browse source 
regtests for CDFT based QMMM simulations (Nico Holmberg)


svn-origin-rev: 17758
This commit is contained in:
Iain Bethune 2017-02-13 12:04:50 +00:00
parent 5624df93cb
commit c869cb8d1b
69 changed files with 5139 additions and 1415 deletions
Download patch file Download diff file Expand all files Collapse all files

40
src/cp_control_types.F
View file

@ -296,7 +296,7 @@ MODULE cp_control_types
should_skip, print_cavity, &
external_control, first_iteration, &
fragment_density, fragments_integrated, &
use_bohr
use_bohr, transfer_pot
REAL(KIND=dp), POINTER, DIMENSION(:) :: coeff, cutoffs, cutoffs_tmp, &
charges_fragment, radii_tmp, &
radii
@ -347,6 +347,7 @@ MODULE cp_control_types
!> \param purge_offset purging is only allowed when more than purge_offset steps have passed since
!> last purge
!> \param istep a counter to keep track of how many steps have passed since the last purge
!> \param ienergy a counter tracking the total number of CDFT energy evaluations
!> \param constraint_type the constraint type to use
!> \param combined_type for combined density+spin constraint, determines which atoms to apply the spin
!> constraint to
@ -366,20 +367,23 @@ MODULE cp_control_types
!> \param reuse_precond logical which determines if a preconditioner can be reused
!> \param purge_history logical which determines if the wfn/constraint history can be purged
!> \param should_purge logical which determines if purging should take place after this CDFT SCF loop
!> \param calculate_metric logical which determines is the ET coupling reliablity metric is computed
!> \param calculate_metric logical which determines if the ET coupling reliablity metric is computed
!> \param mo_coeff save the MO coeffs (for do_et)
!> \param wmat matrix representation of the weight (for do_et)
!> \param wmat matrix representation of the weight function (for do_et)
!> \param matrix_s save the overlap matrix (for do_et)
!> \param matrix_s save the density matrix (for calculate_metric)
!> \param matrix_p save the density matrix (for calculate_metric)
! **************************************************************************************************
! Copied from qs_scf_types to avoid circular dependency
TYPE qs_outer_scf_type
INTEGER :: iter_count
LOGICAL :: deallocate_jacobian
REAL(KIND=dp), DIMENSION(:), POINTER :: energy
REAL(KIND=dp), DIMENSION(:, :), &
POINTER :: variables
REAL(KIND=dp), DIMENSION(:, :), &
POINTER :: gradient
REAL(KIND=dp), DIMENSION(:, :), &
POINTER :: inv_jacobian
INTEGER, DIMENSION(:), POINTER :: count
END TYPE qs_outer_scf_type
@ -387,22 +391,28 @@ MODULE cp_control_types
TYPE outer_scf_control_type
LOGICAL :: have_scf
LOGICAL :: build_jacobian
INTEGER :: max_scf
LOGICAL :: broyden_update
INTEGER :: max_scf, &
ijacobian
REAL(KIND=dp) :: eps_scf, step_size, &
jacobian_step
jacobian_step, &
newton_step
INTEGER :: TYPE
INTEGER :: optimizer
INTEGER :: diis_buffer_length
INTEGER :: extrapolation_order
INTEGER :: bisect_trust_count
INTEGER :: jacobian_type
INTEGER :: broyden_type
INTEGER :: jacobian_freq(2)
END TYPE outer_scf_control_type
TYPE cdft_control_type
INTEGER :: ref_count, total_steps, TYPE, &
precond_freq, nreused, max_reuse, &
purge_freq, nbad_conv, purge_offset, &
istep, constraint_type, combined_type
istep, constraint_type, combined_type, &
ienergy
TYPE(qs_outer_scf_type) :: constraint
TYPE(outer_scf_control_type) :: constraint_control, ot_control
TYPE(hirshfeld_constraint_type), &
@ -889,6 +899,7 @@ CONTAINS
becke_control%natoms = 0
becke_control%need_pot = .TRUE.
becke_control%save_pot = .FALSE.
becke_control%transfer_pot = .FALSE.
becke_control%adjust = .FALSE.
becke_control%atomic_charges = .FALSE.
becke_control%confine_method = becke_none_conf
@ -914,6 +925,9 @@ CONTAINS
becke_control%external_control = .FALSE.
becke_control%constraint_type = cdft_density_constraint
becke_control%combined_type = cdft_combined_all
NULLIFY (becke_control%becke_pot%pw)
NULLIFY (becke_control%cavity%pw)
NULLIFY (becke_control%combined_weight%pw)
NULLIFY (becke_control%strength)
NULLIFY (becke_control%target)
NULLIFY (becke_control%becke_order_p)
@ -1043,6 +1057,7 @@ CONTAINS
cdft_control%istep = 0
cdft_control%constraint_type = cdft_density_constraint
cdft_control%combined_type = cdft_combined_all
cdft_control%ienergy = 0
NULLIFY (cdft_control%hirshfeld_control)
NULLIFY (cdft_control%weight%pw)
NULLIFY (cdft_control%wmat%matrix)
@ -1059,6 +1074,7 @@ CONTAINS
cdft_control%ot_control%diis_buffer_length = -1
cdft_control%ot_control%jacobian_type = -1
cdft_control%ot_control%jacobian_step = 0.0_dp
cdft_control%ot_control%newton_step = 1.0_dp
cdft_control%constraint_control%have_scf = .FALSE.
cdft_control%constraint_control%max_scf = 0
cdft_control%constraint_control%eps_scf = 0.0_dp
@ -1067,12 +1083,20 @@ CONTAINS
cdft_control%constraint_control%optimizer = -1
cdft_control%constraint_control%diis_buffer_length = -1
cdft_control%constraint_control%jacobian_type = -1
cdft_control%constraint_control%broyden_type = -1
cdft_control%constraint_control%jacobian_freq = 1
cdft_control%constraint_control%jacobian_step = 0.0_dp
cdft_control%constraint_control%newton_step = 1.0_dp
cdft_control%constraint_control%ijacobian = 0
cdft_control%constraint_control%build_jacobian = .FALSE.
cdft_control%constraint_control%broyden_update = .FALSE.
cdft_control%constraint%iter_count = 0
NULLIFY (cdft_control%constraint%variables)
NULLIFY (cdft_control%constraint%gradient)
NULLIFY (cdft_control%constraint%energy)
NULLIFY (cdft_control%constraint%count)
NULLIFY (cdft_control%constraint%inv_jacobian)
cdft_control%constraint%deallocate_jacobian = .TRUE.
END SUBROUTINE cdft_control_create
! *****************************************************************************
@ -1116,6 +1140,8 @@ CONTAINS
DEALLOCATE (cdft_control%constraint%gradient)
IF (ASSOCIATED(cdft_control%constraint%energy)) &
DEALLOCATE (cdft_control%constraint%energy)
IF (ASSOCIATED(cdft_control%constraint%inv_jacobian)) &
DEALLOCATE (cdft_control%constraint%inv_jacobian)
DEALLOCATE (cdft_control)
END IF
END SUBROUTINE cdft_control_release

21
src/cp_control_utils.F
View file

@ -1941,6 +1941,7 @@ CONTAINS
CHARACTER(len=*), PARAMETER :: routineN = 'read_cdft_control_section', &
routineP = moduleN//':'//routineN
INTEGER, DIMENSION(:), POINTER :: tmplist
LOGICAL :: exists
TYPE(section_vals_type), POINTER :: hirshfeld_constraint_section, &
outer_scf_section
@ -1961,6 +1962,8 @@ CONTAINS
i_val=qs_control%cdft_control%purge_freq)
CALL section_vals_val_get(cdft_control_section, "PURGE_OFFSET", &
i_val=qs_control%cdft_control%purge_offset)
CALL section_vals_val_get(cdft_control_section, "COUNTER", &
i_val=qs_control%cdft_control%ienergy)
outer_scf_section => section_vals_get_subs_vals(cdft_control_section, "OUTER_SCF")
CALL section_vals_val_get(outer_scf_section, "_SECTION_PARAMETERS_", &
l_val=qs_control%cdft_control%constraint_control%have_scf)
@ -1990,6 +1993,24 @@ CONTAINS
i_val=qs_control%cdft_control%constraint_control%jacobian_type)
CALL section_vals_val_get(outer_scf_section, "JACOBIAN_STEP", &
r_val=qs_control%cdft_control%constraint_control%jacobian_step)
CALL section_vals_val_get(outer_scf_section, "BROYDEN_TYPE", &
i_val=qs_control%cdft_control%constraint_control%broyden_type)
CALL section_vals_val_get(outer_scf_section, "JACOBIAN_FREQ", explicit=exists)
IF (exists) THEN
CALL section_vals_val_get(outer_scf_section, "JACOBIAN_FREQ", &
i_vals=tmplist)
IF (SIZE(tmplist) /= 2) &
CALL cp_abort(__LOCATION__, &
"Keyword JACOBIAN_FREQ takes exactly two input values.")
IF (ANY(tmplist .LT. 0)) &
CALL cp_abort(__LOCATION__, &
"Keyword JACOBIAN_FREQ takes only positive values.")
IF (ALL(tmplist .EQ. 0)) &
CALL cp_abort(__LOCATION__, &
"Both values to keyword JACOBIAN_FREQ cannot be zero.")
qs_control%cdft_control%constraint_control%jacobian_freq(:) = tmplist(1:2)
END IF
SELECT CASE (qs_control%cdft_control%type)
CASE (outer_scf_hirshfeld_constraint)
IF (qs_control%cdft_control%constraint_control%type == outer_scf_cdft_constraint) THEN

3
src/et_coupling.F
View file

@ -1737,7 +1737,8 @@ CONTAINS
END DO
END IF
NULLIFY (gradients)
DEALLOCATE (dft_control%qs_control%becke_control%gradients)
IF (.NOT. dft_control%qs_control%becke_control%transfer_pot) &
DEALLOCATE (dft_control%qs_control%becke_control%gradients)
IF (dft_control%qs_control%becke_control%should_skip .AND. .NOT. &
dft_control%qs_control%becke_control%cavity_confine) &
DEALLOCATE (dft_control%qs_control%becke_control%skip_list)

323
src/force_env_methods.F
View file

@ -82,8 +82,7 @@ MODULE force_env_methods
mix_restrained,&
use_bazant_eip,&
use_lenosky_eip
USE input_section_types, ONLY: section_vals_get,&
section_vals_get_subs_vals,&
USE input_section_types, ONLY: section_vals_get_subs_vals,&
section_vals_retain,&
section_vals_type,&
section_vals_val_get
@ -775,11 +774,11 @@ CONTAINS
CHARACTER(LEN=default_path_length) :: coupling_function
CHARACTER(LEN=default_string_length) :: def_error, description, this_error
INTEGER :: et_freq, iforce_eval, iparticle, &
jparticle, mixing_type, my_group, &
natom, nforce_eval, source, unit_nr
INTEGER :: iforce_eval, iparticle, jparticle, &
mixing_type, my_group, natom, &
nforce_eval, source, unit_nr
INTEGER, DIMENSION(:), POINTER :: glob_natoms, map_index
LOGICAL :: dip_exists, do_mixed_cdft, explicit
LOGICAL :: dip_exists
REAL(KIND=dp) :: coupling_parameter, dedf, der_1, der_2, &
dx, energy, err, lambda, lerr, &
restraint_strength, restraint_target, &
@ -792,12 +791,10 @@ CONTAINS
TYPE(cp_result_type), POINTER :: loc_results, results_mix
TYPE(cp_subsys_p_type), DIMENSION(:), POINTER :: subsystems
TYPE(cp_subsys_type), POINTER :: subsys_mix
TYPE(dft_control_type), POINTER :: dft_control
TYPE(mixed_energy_type), POINTER :: mixed_energy
TYPE(mixed_force_type), DIMENSION(:), POINTER :: global_forces
TYPE(particle_list_p_type), DIMENSION(:), POINTER :: particles
TYPE(particle_list_type), POINTER :: particles_mix
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(section_vals_type), POINTER :: force_env_section, gen_section, &
mapping_section, mixed_section, &
root_section
@ -816,7 +813,7 @@ CONTAINS
particles=particles_mix, &
virial=virial_mix, &
results=results_mix)
NULLIFY (map_index, glob_natoms, global_forces, qs_env, dft_control)
NULLIFY (map_index, glob_natoms, global_forces)
nforce_eval = SIZE(force_env%sub_force_env)
mixed_section => section_vals_get_subs_vals(force_env_section, "MIXED")
@ -833,32 +830,8 @@ CONTAINS
energies = 0.0_dp
glob_natoms = 0
! Check if mixed CDFT calculation is requested and initialize
CALL section_vals_val_get(mixed_section, "MIXING_TYPE", i_val=mixing_type)
IF (mixing_type == mix_linear_combination .AND. &
.NOT. ASSOCIATED(force_env%mixed_env%cdft_control)) THEN
CALL section_vals_val_get(mixed_section, "LINEAR%MIXED_CDFT", l_val=do_mixed_cdft)
force_env%mixed_env%do_mixed_cdft = do_mixed_cdft
IF (force_env%mixed_env%do_mixed_cdft) THEN
! Enforce some limitations
CPASSERT(force_env%mixed_env%ngroups == 2)
CPASSERT(nforce_eval == 2)
CALL section_vals_get(mapping_section, explicit=explicit)
! The sub_force_envs must share the same geometrical structure
CPASSERT(.NOT. explicit)
CALL section_vals_val_get(mixed_section, "LINEAR%MIXED_CDFT_COUPLING", i_val=et_freq)
IF (et_freq .LT. 0) THEN
force_env%mixed_env%do_mixed_et = .FALSE.
ELSE
force_env%mixed_env%do_mixed_et = .TRUE.
IF (et_freq == 0) THEN
force_env%mixed_env%et_freq = 1
ELSE
force_env%mixed_env%et_freq = et_freq
END IF
END IF
CALL mixed_cdft_init(force_env, calculate_forces)
END IF
END IF
CALL mixed_cdft_init(force_env, calculate_forces)
!
IF (.NOT. force_env%mixed_env%do_mixed_cdft) THEN
DO iforce_eval = 1, nforce_eval
NULLIFY (subsystems(iforce_eval)%subsys, particles(iforce_eval)%list)
@ -919,102 +892,13 @@ CONTAINS
CALL cp_rm_default_logger()
END DO
ELSE
DO iforce_eval = 1, nforce_eval
NULLIFY (subsystems(iforce_eval)%subsys, particles(iforce_eval)%list)
NULLIFY (results(iforce_eval)%results, virials(iforce_eval)%virial)
CALL virial_create(virials(iforce_eval)%virial)
CALL cp_result_create(results(iforce_eval)%results)
IF (.NOT. ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env)) CYCLE
! From this point on the error is the sub_error
my_group = force_env%mixed_env%group_distribution(force_env%para_env%mepos)
my_logger => force_env%mixed_env%sub_logger(my_group+1)%p
! Copy iterations info (they are updated only in the main mixed_env)
CALL cp_iteration_info_copy_iter(logger%iter_info, my_logger%iter_info)
! Get all available subsys
CALL force_env_get(force_env=force_env%sub_force_env(iforce_eval)%force_env, &
subsys=subsystems(iforce_eval)%subsys)
! all force_env share the same cell
CALL cp_subsys_set(subsystems(iforce_eval)%subsys, cell=cell_mix)
! Get available particles
CALL cp_subsys_get(subsys=subsystems(iforce_eval)%subsys, &
particles=particles(iforce_eval)%list)
! Get Mapping index array
natom = SIZE(particles(iforce_eval)%list%els)
CALL get_subsys_map_index(mapping_section, natom, iforce_eval, nforce_eval, &
map_index)
! Mapping particles from iforce_eval environment to the mixed env
DO iparticle = 1, natom
jparticle = map_index(iparticle)
particles(iforce_eval)%list%els(iparticle)%r = particles_mix%els(jparticle)%r
END DO
! Mixed CDFT + QMMM: Need to translate now
IF (force_env%mixed_env%do_mixed_qmmm_cdft) &
CALL apply_qmmm_translate(force_env%sub_force_env(iforce_eval)%force_env%qmmm_env)
END DO
! For mixed CDFT, build weight and gradient on all processors before splitting into groups and
! starting energy calculation
IF (force_env%mixed_env%do_mixed_cdft) THEN
CALL mixed_cdft_build_weight(force_env, calculate_forces)
END IF
DO iforce_eval = 1, nforce_eval
IF (.NOT. ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env)) CYCLE
CALL cp_add_default_logger(my_logger)
! Calculate energy and forces for each sub_force_env
CALL force_env_calc_energy_force(force_env%sub_force_env(iforce_eval)%force_env, &
calc_force=calculate_forces, &
skip_external_control=.TRUE.)
! Only the rank 0 process collect info for each computation
IF (force_env%sub_force_env(iforce_eval)%force_env%para_env%mepos == &
force_env%sub_force_env(iforce_eval)%force_env%para_env%source) THEN
CALL force_env_get(force_env%sub_force_env(iforce_eval)%force_env, &
potential_energy=energy)
CALL cp_subsys_get(subsystems(iforce_eval)%subsys, &
virial=loc_virial, results=loc_results)
energies(iforce_eval) = energy
glob_natoms(iforce_eval) = natom
CALL cp_virial(loc_virial, virials(iforce_eval)%virial)
CALL cp_result_copy(loc_results, results(iforce_eval)%results)
END IF
! Deallocate map_index array
IF (ASSOCIATED(map_index)) THEN
DEALLOCATE (map_index)
END IF
CALL cp_rm_default_logger()
END DO
CALL mixed_cdft_energy_forces(force_env, calculate_forces, particles, energies, &
glob_natoms, virials, results)
END IF
! Handling Parallel execution
CALL mp_sync(force_env%para_env%group)
! Transfer cdft strengths for writing restart
IF (force_env%mixed_env%do_mixed_cdft) THEN
IF (.NOT. ASSOCIATED(force_env%mixed_env%strength)) &
ALLOCATE (force_env%mixed_env%strength(2))
force_env%mixed_env%strength = 0.0_dp
DO iforce_eval = 1, nforce_eval
IF (.NOT. ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env)) CYCLE
IF (force_env%mixed_env%do_mixed_qmmm_cdft) THEN
qs_env => force_env%sub_force_env(iforce_eval)%force_env%qmmm_env%qs_env
ELSE
CALL force_env_get(force_env%sub_force_env(iforce_eval)%force_env, qs_env=qs_env)
END IF
CALL get_qs_env(qs_env, dft_control=dft_control)
IF (force_env%sub_force_env(iforce_eval)%force_env%para_env%mepos == &
force_env%sub_force_env(iforce_eval)%force_env%para_env%source) &
force_env%mixed_env%strength(iforce_eval) = dft_control%qs_control%cdft_control%strength(1)
END DO
CALL mp_sum(force_env%mixed_env%strength, force_env%para_env%group)
END IF
! Mixed CDFT: calculate ET coupling
IF (force_env%mixed_env%do_mixed_et) THEN
IF (MODULO(force_env%mixed_env%cdft_control%sim_step, force_env%mixed_env%et_freq) == 0) &
CALL mixed_cdft_calculate_coupling(force_env)
END IF
! Post CDFT operations
CALL mixed_cdft_post_energy_forces(force_env)
! Let's transfer energy, natom, forces, virials
CALL mp_sum(energies, force_env%para_env%group)
CALL mp_sum(glob_natoms, force_env%para_env%group)
@ -1184,4 +1068,187 @@ CONTAINS
CALL cp_print_key_finished_output(unit_nr, logger, mixed_section, "PRINT%DIPOLE")
END SUBROUTINE mixed_energy_forces
! **************************************************************************************************
!> \brief Driver routine for mixed CDFT energy and force calculations
!> \param force_env the force_env that holds the mixed_env
!> \param calculate_forces if forces should be calculated
!> \param particles system particles
!> \param energies the energies of the CDFT states
!> \param glob_natoms the total number of particles
!> \param virials the virials stored in subsys
!> \param results results stored in subsys
!> \par History
!> 01.17 created [Nico Holmberg]
!> \author Nico Holmberg
! **************************************************************************************************
SUBROUTINE mixed_cdft_energy_forces(force_env, calculate_forces, particles, energies, &
glob_natoms, virials, results)
TYPE(force_env_type), POINTER :: force_env
LOGICAL, INTENT(IN) :: calculate_forces
TYPE(particle_list_p_type), DIMENSION(:), POINTER :: particles
REAL(KIND=dp), DIMENSION(:), POINTER :: energies
INTEGER, DIMENSION(:), POINTER :: glob_natoms
TYPE(virial_p_type), DIMENSION(:), POINTER :: virials
TYPE(cp_result_p_type), DIMENSION(:), POINTER :: results
CHARACTER(len=*), PARAMETER :: routineN = 'mixed_cdft_energy_forces', &
routineP = moduleN//':'//routineN
INTEGER :: iforce_eval, iparticle, jparticle, &
my_group, natom, nforce_eval
INTEGER, DIMENSION(:), POINTER :: map_index
REAL(KIND=dp) :: energy
TYPE(cell_type), POINTER :: cell_mix
TYPE(cp_logger_type), POINTER :: logger, my_logger
TYPE(cp_result_type), POINTER :: loc_results, results_mix
TYPE(cp_subsys_p_type), DIMENSION(:), POINTER :: subsystems
TYPE(cp_subsys_type), POINTER :: subsys_mix
TYPE(particle_list_type), POINTER :: particles_mix
TYPE(section_vals_type), POINTER :: force_env_section, mapping_section, &
mixed_section, root_section
TYPE(virial_type), POINTER :: loc_virial, virial_mix
logger => cp_get_default_logger()
CPASSERT(ASSOCIATED(force_env))
! Get infos about the mixed subsys
CALL force_env_get(force_env=force_env, &
subsys=subsys_mix, &
force_env_section=force_env_section, &
root_section=root_section, &
cell=cell_mix)
CALL cp_subsys_get(subsys=subsys_mix, &
particles=particles_mix, &
virial=virial_mix, &
results=results_mix)
NULLIFY (map_index)
nforce_eval = SIZE(force_env%sub_force_env)
mixed_section => section_vals_get_subs_vals(force_env_section, "MIXED")
mapping_section => section_vals_get_subs_vals(mixed_section, "MAPPING")
ALLOCATE (subsystems(nforce_eval))
DO iforce_eval = 1, nforce_eval
NULLIFY (subsystems(iforce_eval)%subsys, particles(iforce_eval)%list)
NULLIFY (results(iforce_eval)%results, virials(iforce_eval)%virial)
CALL virial_create(virials(iforce_eval)%virial)
CALL cp_result_create(results(iforce_eval)%results)
IF (.NOT. ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env)) CYCLE
! Get all available subsys
CALL force_env_get(force_env=force_env%sub_force_env(iforce_eval)%force_env, &
subsys=subsystems(iforce_eval)%subsys)
! all force_env share the same cell
CALL cp_subsys_set(subsystems(iforce_eval)%subsys, cell=cell_mix)
! Get available particles
CALL cp_subsys_get(subsys=subsystems(iforce_eval)%subsys, &
particles=particles(iforce_eval)%list)
! Get Mapping index array
natom = SIZE(particles(iforce_eval)%list%els)
! Serial mode need to deallocate first
IF (ASSOCIATED(map_index)) &
DEALLOCATE (map_index)
CALL get_subsys_map_index(mapping_section, natom, iforce_eval, nforce_eval, &
map_index)
! Mapping particles from iforce_eval environment to the mixed env
DO iparticle = 1, natom
jparticle = map_index(iparticle)
particles(iforce_eval)%list%els(iparticle)%r = particles_mix%els(jparticle)%r
END DO
! Mixed CDFT + QMMM: Need to translate now
IF (force_env%mixed_env%do_mixed_qmmm_cdft) &
CALL apply_qmmm_translate(force_env%sub_force_env(iforce_eval)%force_env%qmmm_env)
END DO
! For mixed CDFT calculations parallized over CDFT states
! build weight and gradient on all processors before splitting into groups and
! starting energy calculation
CALL mixed_cdft_build_weight(force_env, calculate_forces)
DO iforce_eval = 1, nforce_eval
IF (.NOT. ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env)) CYCLE
! From this point on the error is the sub_error
IF (.NOT. force_env%mixed_env%cdft_control%is_parallel .AND. iforce_eval == 2) THEN
my_logger => force_env%mixed_env%cdft_control%sub_logger
ELSE
my_group = force_env%mixed_env%group_distribution(force_env%para_env%mepos)
my_logger => force_env%mixed_env%sub_logger(my_group+1)%p
END IF
! Copy iterations info (they are updated only in the main mixed_env)
CALL cp_iteration_info_copy_iter(logger%iter_info, my_logger%iter_info)
CALL cp_add_default_logger(my_logger)
! Serial CDFT calculation: transfer weight/gradient
CALL mixed_cdft_build_weight(force_env, calculate_forces, iforce_eval)
! Calculate energy and forces for each sub_force_env
CALL force_env_calc_energy_force(force_env%sub_force_env(iforce_eval)%force_env, &
calc_force=calculate_forces, &
skip_external_control=.TRUE.)
! Only the rank 0 process collect info for each computation
IF (force_env%sub_force_env(iforce_eval)%force_env%para_env%mepos == &
force_env%sub_force_env(iforce_eval)%force_env%para_env%source) THEN
CALL force_env_get(force_env%sub_force_env(iforce_eval)%force_env, &
potential_energy=energy)
CALL cp_subsys_get(subsystems(iforce_eval)%subsys, &
virial=loc_virial, results=loc_results)
energies(iforce_eval) = energy
glob_natoms(iforce_eval) = natom
CALL cp_virial(loc_virial, virials(iforce_eval)%virial)
CALL cp_result_copy(loc_results, results(iforce_eval)%results)
END IF
! Deallocate map_index array
IF (ASSOCIATED(map_index)) THEN
DEALLOCATE (map_index)
END IF
CALL cp_rm_default_logger()
END DO
DEALLOCATE (subsystems)
END SUBROUTINE mixed_cdft_energy_forces
! **************************************************************************************************
!> \brief Perform additional tasks for mixed CDFT calculations after solving the electronic structure
!> of both CDFT states
!> \param force_env the force_env that holds the CDFT states
!> \par History
!> 01.17 created [Nico Holmberg]
!> \author Nico Holmberg
! **************************************************************************************************
SUBROUTINE mixed_cdft_post_energy_forces(force_env)
TYPE(force_env_type), POINTER :: force_env
CHARACTER(len=*), PARAMETER :: routineN = 'mixed_cdft_post_energy_forces', &
routineP = moduleN//':'//routineN
INTEGER :: iforce_eval, nforce_eval
TYPE(dft_control_type), POINTER :: dft_control
TYPE(qs_environment_type), POINTER :: qs_env
CPASSERT(ASSOCIATED(force_env))
NULLIFY (qs_env, dft_control)
IF (force_env%mixed_env%do_mixed_cdft) THEN
nforce_eval = SIZE(force_env%sub_force_env)
! Transfer cdft strengths for writing restart
IF (.NOT. ASSOCIATED(force_env%mixed_env%strength)) &
ALLOCATE (force_env%mixed_env%strength(2))
force_env%mixed_env%strength = 0.0_dp
DO iforce_eval = 1, nforce_eval
IF (.NOT. ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env)) CYCLE
IF (force_env%mixed_env%do_mixed_qmmm_cdft) THEN
qs_env => force_env%sub_force_env(iforce_eval)%force_env%qmmm_env%qs_env
ELSE
CALL force_env_get(force_env%sub_force_env(iforce_eval)%force_env, qs_env=qs_env)
END IF
CALL get_qs_env(qs_env, dft_control=dft_control)
IF (force_env%sub_force_env(iforce_eval)%force_env%para_env%mepos == &
force_env%sub_force_env(iforce_eval)%force_env%para_env%source) &
force_env%mixed_env%strength(iforce_eval) = dft_control%qs_control%cdft_control%strength(1)
END DO
CALL mp_sum(force_env%mixed_env%strength, force_env%para_env%group)
! Mixed CDFT: calculate ET coupling
IF (force_env%mixed_env%do_mixed_et) THEN
IF (MODULO(force_env%mixed_env%cdft_control%sim_step, force_env%mixed_env%et_freq) == 0) &
CALL mixed_cdft_calculate_coupling(force_env)
END IF
END IF
END SUBROUTINE mixed_cdft_post_energy_forces
END MODULE force_env_methods

14
src/input_constants.F
View file

@ -647,13 +647,19 @@ MODULE input_constants
outer_scf_optimizer_none = 1003, &
outer_scf_optimizer_bisect = 1004, &
outer_scf_optimizer_broyden = 1005, &
outer_scf_optimizer_newton = 1006
outer_scf_optimizer_newton = 1006, &
outer_scf_optimizer_secant = 1007, &
outer_scf_optimizer_newton_ls = 1008
! outer scf broyden optimizer types
INTEGER, PARAMETER, PUBLIC :: broyden_type_1 = 1101, &
broyden_type_2 = 1102, &
broyden_type_3 = 1103, &
broyden_type_4 = 1104
broyden_type_1_explicit = 1102, &
broyden_type_2 = 1103, &
broyden_type_2_explicit = 1104, &
broyden_type_1_ls = 1105, &
broyden_type_1_explicit_ls = 1106, &
broyden_type_2_ls = 1107, &
broyden_type_2_explicit_ls = 1108
! finite difference types for calculation of inverse jacobian
INTEGER, PARAMETER, PUBLIC :: jacobian_fd1 = 1, &

107
src/input_cp2k_dft.F
View file

@ -30,7 +30,9 @@ MODULE input_cp2k_dft
cp_unit_to_cp2k
USE input_constants, ONLY: &
atomic_guess, becke_cavity_conf, becke_cutoff_element, becke_cutoff_global, &
becke_dynamic_conf, becke_mixed_conf, becke_none_conf, becke_static_conf, casci_canonical, &
becke_dynamic_conf, becke_mixed_conf, becke_none_conf, becke_static_conf, broyden_type_1, &
broyden_type_1_explicit, broyden_type_1_explicit_ls, broyden_type_1_ls, broyden_type_2, &
broyden_type_2_explicit, broyden_type_2_explicit_ls, broyden_type_2_ls, casci_canonical, &
cdft_combined_acceptor, cdft_combined_all, cdft_combined_constraint, cdft_combined_donor, &
cdft_density_constraint, cdft_magnetization_constraint, cholesky_dbcsr, cholesky_inverse, &
cholesky_off, cholesky_reduce, cholesky_restore, core_guess, diag_block_davidson, &
@ -70,24 +72,25 @@ MODULE input_cp2k_dft
outer_scf_becke_constraint, outer_scf_cdft_constraint, outer_scf_ddapc_constraint, &
outer_scf_hirshfeld_constraint, outer_scf_none, outer_scf_optimizer_bisect, &
outer_scf_optimizer_broyden, outer_scf_optimizer_diis, outer_scf_optimizer_newton, &
outer_scf_optimizer_none, outer_scf_optimizer_sd, outer_scf_s2_constraint, plus_u_lowdin, &
plus_u_mulliken, plus_u_mulliken_charges, pw_interp, radius_covalent, radius_default, &
radius_single, radius_user, radius_vdw, random_guess, real_time_propagation, &
ref_charge_atomic, ref_charge_mulliken, rel_dkh, rel_none, rel_pot_erfc, rel_pot_full, &
rel_sczora_mp, rel_trans_atom, rel_trans_full, rel_trans_molecule, rel_zora, &
rel_zora_full, rel_zora_mp, restart_guess, rsdft_model, sccs_andreussi, &
sccs_derivative_cd3, sccs_derivative_cd5, sccs_derivative_cd7, sccs_derivative_fft, &
sccs_fattebert_gygi, shape_function_density, shape_function_gaussian, sic_ad, sic_eo, &
sic_list_all, sic_list_unpaired, sic_mauri_spz, sic_mauri_us, sic_none, slater, &
smear_energy_window, smear_fermi_dirac, smear_list, sparse_guess, spline3_nopbc_interp, &
spline3_pbc_interp, tddfpt_davidson, tddfpt_excitations, tddfpt_lanczos, tddfpt_singlet, &
tddfpt_triplet, use_coulomb, use_diff, use_no, use_restart_wfn, use_rt_restart, &
use_scf_wfn, wannier_projection, weight_type_mass, weight_type_unit, wfi_aspc_nr, &
wfi_frozen_method_nr, wfi_linear_p_method_nr, wfi_linear_ps_method_nr, &
wfi_linear_wf_method_nr, wfi_ps_method_nr, wfi_use_guess_method_nr, &
wfi_use_prev_p_method_nr, wfi_use_prev_rho_r_method_nr, wfi_use_prev_wf_method_nr, &
xas_1s_type, xas_2p_type, xas_2s_type, xas_dip_len, xas_dip_vel, xas_dscf, xas_none, &
xas_tp_fh, xas_tp_flex, xas_tp_hh, xas_tp_xfh, xas_tp_xhh, xes_tp_val
outer_scf_optimizer_newton_ls, outer_scf_optimizer_none, outer_scf_optimizer_sd, &
outer_scf_optimizer_secant, outer_scf_s2_constraint, plus_u_lowdin, plus_u_mulliken, &
plus_u_mulliken_charges, pw_interp, radius_covalent, radius_default, radius_single, &
radius_user, radius_vdw, random_guess, real_time_propagation, ref_charge_atomic, &
ref_charge_mulliken, rel_dkh, rel_none, rel_pot_erfc, rel_pot_full, rel_sczora_mp, &
rel_trans_atom, rel_trans_full, rel_trans_molecule, rel_zora, rel_zora_full, rel_zora_mp, &
restart_guess, rsdft_model, sccs_andreussi, sccs_derivative_cd3, sccs_derivative_cd5, &
sccs_derivative_cd7, sccs_derivative_fft, sccs_fattebert_gygi, shape_function_density, &
shape_function_gaussian, sic_ad, sic_eo, sic_list_all, sic_list_unpaired, sic_mauri_spz, &
sic_mauri_us, sic_none, slater, smear_energy_window, smear_fermi_dirac, smear_list, &
sparse_guess, spline3_nopbc_interp, spline3_pbc_interp, tddfpt_davidson, &
tddfpt_excitations, tddfpt_lanczos, tddfpt_singlet, tddfpt_triplet, use_coulomb, use_diff, &
use_no, use_restart_wfn, use_rt_restart, use_scf_wfn, wannier_projection, &
weight_type_mass, weight_type_unit, wfi_aspc_nr, wfi_frozen_method_nr, &
wfi_linear_p_method_nr, wfi_linear_ps_method_nr, wfi_linear_wf_method_nr, &
wfi_ps_method_nr, wfi_use_guess_method_nr, wfi_use_prev_p_method_nr, &
wfi_use_prev_rho_r_method_nr, wfi_use_prev_wf_method_nr, xas_1s_type, xas_2p_type, &
xas_2s_type, xas_dip_len, xas_dip_vel, xas_dscf, xas_none, xas_tp_fh, xas_tp_flex, &
xas_tp_hh, xas_tp_xfh, xas_tp_xhh, xes_tp_val
USE input_cp2k_almo, ONLY: create_almo_scf_section
USE input_cp2k_distribution, ONLY: create_distribution_section
USE input_cp2k_field, ONLY: create_efield_section,&
@ -4340,7 +4343,7 @@ CONTAINS
description="Parameters needed to set up a CDFT calculation with OT."// &
" Constraint is converged in a separate external SCF loop with settings"// &
" read from the OUTER_SCF section. Supported constraints: Gaussian Hirshfeld (partial)"// &
" Becke.", n_keywords=7, n_subsections=2, &
" Becke.", n_keywords=8, n_subsections=2, &
repeats=.FALSE., citations=(/Holmberg2017/))
NULLIFY (subsection, keyword)
@ -4408,6 +4411,13 @@ CONTAINS
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, name="COUNTER", &
description="A counter to track the total number of energy evaluations. Needed by"// &
" some optimizers to print information. Useful for restarts.", &
usage="COUNTER {int} ", default_i_val=0)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
END SUBROUTINE create_cdft_control_section
! *****************************************************************************
@ -5339,7 +5349,7 @@ CONTAINS
CPASSERT(.NOT. ASSOCIATED(section))
CALL section_create(section, "OUTER_SCF", &
description="parameters controlling the outer SCF loop", &
n_keywords=11, n_subsections=0, repeats=.FALSE.)
n_keywords=13, n_subsections=0, repeats=.FALSE.)
NULLIFY (keyword)
@ -5359,12 +5369,12 @@ CONTAINS
"Enforce a constraint on the S2, requires the corresponding section", &
"Enforce a constraint on the Becke weight population, "// &
"requires the corresponding section", &
"Optimize positions of basis functions, if atom types FLOATING_BASIS_CENTER "// &
" are defined", &
"Enforce a constraint on a generic CDFT weight population, "// &
"requires the corresponding section"// &
" which determines the type of weight used "// &
"(currently only Gaussian Hirshfeld supported)", &
"Optimize positions of basis functions, if atom types FLOATING_BASIS_CENTER "// &
" are defined", &
"Do nothing in the outer loop, useful for resetting the inner loop,"), &
enum_i_vals=(/outer_scf_ddapc_constraint, outer_scf_s2_constraint, &
outer_scf_becke_constraint, outer_scf_basis_center_opt, &
@ -5376,22 +5386,54 @@ CONTAINS
description="Method used to bring the outer loop to a stationary point", &
usage="OPTIMIZER SD", &
default_i_val=outer_scf_optimizer_none, &
enum_c_vals=s2a("SD", "DIIS", "NONE", "BISECT", "BROYDEN", "NEWTON"), &
enum_c_vals=s2a("SD", "DIIS", "NONE", "BISECT", "BROYDEN", "NEWTON", "SECANT", "NEWTON_LS"), &
enum_desc=s2a("Takes steps in the direction of the gradient, multiplied by step_size", &
"Uses a Direct Inversion in the Iterative Subspace method", &
"Do nothing, useful only with the none type", &
"Bisection of the gradient, useful for difficult one dimensional cases", &
"Broyden's method. Variant defined in BROYDEN_TYPE.", &
"Newton's method."), &
"Newton's method. Only compatible with CDFT constraints", &
"Secant method. Only for one dimensional cases. See Broyden for "// &
"multidimensional cases.", &
"Newton's method with backtracking line search to find the optimal step size. "// &
"Only compatible with CDFT constraints. Starts from the regular Newton solution "// &
"and successively reduces the step size until the L2 norm of the CDFT gradient "// &
"decreases. Potentially very expensive because each iteration performs a full "// &
" SCF calculation."), &
enum_i_vals=(/outer_scf_optimizer_sd, outer_scf_optimizer_diis, outer_scf_optimizer_none, &
outer_scf_optimizer_bisect, outer_scf_optimizer_broyden, &
outer_scf_optimizer_newton/))
outer_scf_optimizer_newton, outer_scf_optimizer_secant, &
outer_scf_optimizer_newton_ls/))
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, name="BROYDEN_TYPE", &
description="Specifies the Broyden optimizer variant to use.", &
usage="BROYDEN_TYPE BT1", &
default_i_val=broyden_type_1, &
enum_c_vals=s2a("BT1", "BT1_EXPLICIT", "BT2", "BT2_EXPLICIT", &
"BT1_LS", "BT1_EXPLICIT_LS", "BT2_LS", "BT2_EXPLICIT_LS"), &
enum_desc=s2a("Broyden's first method, also known as the good method. The initial Jacobian"// &
" is built from MD history if available. Otherwise switches to SD for one"// &
" SCF iteration until a Jacobian can be built from the SCF history.", &
"Same as BT1, but computes the explicit Jacobian with finite differences. "// &
"Requires a CDFT SCF procedure to be active.", &
"Same as BT1, but uses Broyden's second method, also known as the bad method", &
"Same as BT1_EXPLICIT, but using Broyden's second method.", &
"Same as BT1, but uses backtracking line search for optimizing the step size "// &
"(see optimizer NEWTON_LS).", &
"Same as BT1_EXPLICIT, but uses backtracking line search for optimizing the step size ", &
"Same as BT2, but uses backtracking line search for optimizing the step size ", &
"Same as BT2_EXPLICIT, but uses backtracking line search for optimizing the step size "), &
enum_i_vals=(/broyden_type_1, broyden_type_1_explicit, broyden_type_2, &
broyden_type_2_explicit, broyden_type_1_ls, broyden_type_1_explicit_ls, &
broyden_type_2_ls, broyden_type_2_explicit_ls/))
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, name="JACOBIAN_TYPE", &
description="Finite difference method used to calculate the inverse Jacobian "// &
"needed by some optimizers.", &
"needed by some optimizers. Compatible only with CDFT constraints.", &
usage="JACOBIAN_TYPE FD1", &
default_i_val=jacobian_fd1, &
enum_c_vals=s2a("FD1", "FD1_BACKWARD", "FD2", "FD2_BACKWARD", "FD1_CENTRAL"), &
@ -5411,6 +5453,17 @@ CONTAINS
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, name="JACOBIAN_FREQ", &
description="Defines parameters that control how often the explicit Jacobian is built, "// &
"which is needed by some optimizers. Expects two values. The first value "// &
"determines how many consecutive CDFT SCF iterations should skip a rebuild,"// &
" whereas the latter how many MD steps. The values can be zero (meaning never "// &
" rebuild) or positive. Both values cannot be zero.", &
usage="JACOBIAN_FREQ 1 1", n_var=2, &
default_i_vals=(/1, 1/), type_of_var=integer_t)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, name="BISECT_TRUST_COUNT", &
description="Maximum number of times the same point will be used in bisection,"// &
" a small number guards against the effect of wrongly converged states.", &

25
src/input_restart_force_eval.F
View file

@ -8,6 +8,7 @@ MODULE input_restart_force_eval
USE atomic_kind_types, ONLY: get_atomic_kind
USE cell_types, ONLY: cell_type,&
real_to_scaled
USE cp_control_types, ONLY: dft_control_type
USE cp_files, ONLY: close_file,&
open_file
USE cp_linked_list_val, ONLY: cp_sll_val_create,&
@ -90,14 +91,16 @@ CONTAINS
INTEGER :: i_subsys, iforce_eval, myid, nforce_eval
INTEGER, DIMENSION(:), POINTER :: i_force_eval
LOGICAL :: multiple_subsys, skip_vel_section
REAL(KIND=dp), DIMENSION(:), POINTER :: work
TYPE(cell_type), POINTER :: cell
TYPE(cp_subsys_type), POINTER :: subsys
TYPE(dft_control_type), POINTER :: dft_control
TYPE(section_vals_type), POINTER :: cell_section, dft_section, &
force_env_sections, qmmm_section, &
rng_section, subsys_section
TYPE(virial_type), POINTER :: virial
NULLIFY (rng_section, subsys_section, cell_section, virial, subsys, cell)
NULLIFY (rng_section, subsys_section, cell_section, virial, subsys, cell, dft_control, work)
! If it's not a dynamical run don't update the velocity section
CALL section_vals_val_get(root_section, "GLOBAL%RUN_TYPE", i_val=myid)
skip_vel_section = ((myid /= mol_dyn_run) .AND. &
@ -136,6 +139,24 @@ CONTAINS
i_rep_section=i_force_eval(iforce_eval))
CALL update_qmmm(qmmm_section, force_env)
! Non-mixed CDFT calculation: update constraint Lagrangian and counter
IF (nforce_eval == 1) THEN
IF (ASSOCIATED(force_env%qs_env)) THEN
CALL get_qs_env(force_env%qs_env, dft_control=dft_control)
IF (dft_control%qs_control%becke_restraint) THEN
ALLOCATE (work(SIZE(dft_control%qs_control%becke_control%strength)))
work = dft_control%qs_control%becke_control%strength
dft_section => section_vals_get_subs_vals3(force_env_sections, "DFT", &
i_rep_section=i_force_eval(iforce_eval))
CALL section_vals_val_set(dft_section, "QS%BECKE_RESTRAINT%STRENGTH", &
r_vals_ptr=work)
END IF
IF (dft_control%qs_control%cdft) THEN
CALL section_vals_val_set(dft_section, "QS%CDFT%COUNTER", &
i_val=dft_control%qs_control%cdft_control%ienergy)
END IF
END IF
END IF
! And now update only the cells of all other force_eval
! This is to make consistent for cell variable runs..
IF (nforce_eval > 1) THEN
@ -158,6 +179,8 @@ CONTAINS
i_rep_section=i_force_eval(iforce_eval))
CALL section_vals_val_set(dft_section, "QS%BECKE_RESTRAINT%STRENGTH", &
r_val=force_env%mixed_env%strength(iforce_eval-1))
CALL section_vals_val_set(dft_section, "QS%CDFT%COUNTER", &
i_val=force_env%mixed_env%cdft_control%sim_step)
END DO
END IF
END IF

1339
src/mixed_cdft_methods.F
View file

File diff suppressed because it is too large Load diff

41
src/mixed_cdft_types.F
View file

@ -14,6 +14,8 @@ MODULE mixed_cdft_types
cp_blacs_env_type
USE cp_control_types, ONLY: becke_control_release,&
becke_restraint_type
USE cp_log_handling, ONLY: cp_logger_release,&
cp_logger_type
USE kinds, ONLY: dp
USE pw_env_types, ONLY: pw_env_release,&
pw_env_type
@ -152,7 +154,8 @@ MODULE mixed_cdft_types
is_special, first_iteration, &
calculate_metric, &
wfn_overlap_method, &
has_unit_metric
has_unit_metric, &
is_parallel
REAL(KIND=dp) :: eps_rho_rspace, sim_dt, &
eps_svd
REAL(KIND=dp), POINTER, DIMENSION(:, :, :) :: weight, cavity
@ -163,11 +166,39 @@ MODULE mixed_cdft_types
TYPE(qs_kind_type), DIMENSION(:), &
POINTER :: qs_kind_set
TYPE(pw_env_type), POINTER :: pw_env
TYPE(cp_logger_type), POINTER :: sub_logger
END TYPE mixed_cdft_type
! **************************************************************************************************
!> \brief Container for constraint settings to check consistency of force_evals
! **************************************************************************************************
TYPE mixed_cdft_settings_type
LOGICAL :: is_spherical, &
is_odd
INTEGER :: nbecke, &
max_nkinds
INTEGER, DIMENSION(2, 3) :: bo
INTEGER, DIMENSION(:), POINTER :: grid_span, &
spherical, &
odd
INTEGER, DIMENSION(:, :), POINTER :: si, &
rs_dims, &
atoms, &
npts
REAL(KIND=dp) :: radius
REAL(KIND=dp), DIMENSION(:), POINTER :: cutoff, &
rel_cutoff
REAL(KIND=dp), DIMENSION(:, :), POINTER :: sr, &
coeffs, &
confine_bounds, &
cutoffs, &
radii
END TYPE mixed_cdft_settings_type
! *** Public data types ***
PUBLIC :: mixed_cdft_type
PUBLIC :: mixed_cdft_type, &
mixed_cdft_settings_type
! *** Public subroutines ***
@ -195,6 +226,7 @@ CONTAINS
NULLIFY (cdft_control%source_list_save, cdft_control%source_bo_save, cdft_control%source_list_bo)
NULLIFY (cdft_control%cavity, cdft_control%weight, cdft_control%sendbuff)
NULLIFY (cdft_control%becke_control, cdft_control%recv_bo)
NULLIFY (cdft_control%sub_logger)
END SUBROUTINE mixed_cdft_type_create
@ -242,11 +274,14 @@ CONTAINS
END DO
DEALLOCATE (cdft_control%sendbuff)
END IF
CALL becke_control_release(cdft_control%becke_control)
IF (ASSOCIATED(cdft_control%becke_control)) &
CALL becke_control_release(cdft_control%becke_control)
IF (ASSOCIATED(cdft_control%blacs_env)) &
CALL cp_blacs_env_release(cdft_control%blacs_env)
IF (ASSOCIATED(cdft_control%qs_kind_set)) &
CALL deallocate_qs_kind_set(cdft_control%qs_kind_set)
IF (ASSOCIATED(cdft_control%sub_logger)) &
CALL cp_logger_release(cdft_control%sub_logger)
DEALLOCATE (cdft_control)
END SUBROUTINE mixed_cdft_type_release

1167
src/mixed_cdft_utils.F Normal file
View file

File diff suppressed because it is too large Load diff

4
src/qs_energy_init.F
View file

@ -251,8 +251,8 @@ CONTAINS
CALL get_qs_env(qs_env, sab_kp=sab_nl)
CALL kpoint_init_cell_index(kpoints, sab_nl, para_env, dft_control)
ENDIF
IF (.NOT. dft_control%qs_control%becke_control%external_control) &
IF (.NOT. (dft_control%qs_control%becke_control%external_control .OR. &
ASSOCIATED(dft_control%qs_control%becke_control%becke_pot%pw))) &
dft_control%qs_control%becke_control%need_pot = .TRUE.
IF (dft_control%qs_control%cdft) THEN
IF (dft_control%qs_control%cdft_control%type == outer_scf_hirshfeld_constraint) &

11
src/qs_environment_methods.F
View file

@ -91,7 +91,8 @@ CONTAINS
INTEGER :: handle, n_ao_aux_fit, nhistory, &
nvariables
REAL(KIND=dp), DIMENSION(:, :), POINTER :: outer_scf_history
REAL(KIND=dp), DIMENSION(:, :), POINTER :: gradient_history, outer_scf_history, &
variable_history
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
TYPE(cell_type), POINTER :: cell
TYPE(cp_blacs_env_type), POINTER :: blacs_env
@ -140,7 +141,13 @@ CONTAINS
nvariables = outer_loop_variables_count(scf_control)
nhistory = scf_control%outer_scf%extrapolation_order
ALLOCATE (outer_scf_history(nvariables, nhistory))
CALL set_qs_env(qs_env, outer_scf_history=outer_scf_history)
ALLOCATE (gradient_history(nvariables, 2))
gradient_history = 0.0_dp
ALLOCATE (variable_history(nvariables, 2))
variable_history = 0.0_dp
CALL set_qs_env(qs_env, outer_scf_history=outer_scf_history, &
gradient_history=gradient_history, &
variable_history=variable_history)
CALL set_qs_env(qs_env, outer_scf_ihistory=0)
ENDIF

23
src/qs_environment_types.F
View file

@ -260,6 +260,8 @@ MODULE qs_environment_types
TYPE(cp_ddapc_ewald_type), POINTER :: cp_ddapc_ewald
REAL(KIND=dp), DIMENSION(:, :), POINTER :: outer_scf_history
INTEGER :: outer_scf_ihistory
REAL(KIND=dp), DIMENSION(:, :), POINTER :: gradient_history, &
variable_history
TYPE(hfx_type), DIMENSION(:, :), POINTER :: x_data
TYPE(et_coupling_type), POINTER :: et_coupling
TYPE(qs_dftb_pairpot_type), DIMENSION(:, :), POINTER :: dftb_potential
@ -456,6 +458,8 @@ CONTAINS
!> \param forces_up_to_date ...
!> \param mscfg_env ...
!> \param almo_scf_env ...
!> \param gradient_history ...
!> \param variable_history ...
!> \date 23.01.2002
!> \author MK
!> \version 1.0
@ -485,7 +489,8 @@ CONTAINS
se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, admm_dm, &
lri_env, lri_density, dispersion_env, vee, rho_external, external_vxc, mask, &
mp2_env, kg_env, WannierCentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, &
s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env)
s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, &
gradient_history, variable_history)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(atomic_kind_type), DIMENSION(:), OPTIONAL, &
POINTER :: atomic_kind_set
@ -613,6 +618,7 @@ CONTAINS
TYPE(molecular_scf_guess_env_type), OPTIONAL, &
POINTER :: mscfg_env
TYPE(almo_scf_env_type), OPTIONAL, POINTER :: almo_scf_env
REAL(KIND=dp), DIMENSION(:, :), OPTIONAL, POINTER :: gradient_history, variable_history
CHARACTER(len=*), PARAMETER :: routineN = 'get_qs_env', routineP = moduleN//':'//routineN
@ -625,6 +631,8 @@ CONTAINS
IF (PRESENT(outer_scf_history)) outer_scf_history => qs_env%outer_scf_history
IF (PRESENT(outer_scf_ihistory)) outer_scf_ihistory = qs_env%outer_scf_ihistory
IF (PRESENT(gradient_history)) gradient_history => qs_env%gradient_history
IF (PRESENT(variable_history)) variable_history => qs_env%variable_history
IF (PRESENT(mp2_env)) mp2_env => qs_env%mp2_env
IF (PRESENT(kg_env)) kg_env => qs_env%kg_env
IF (PRESENT(super_cell)) super_cell => qs_env%super_cell
@ -845,6 +853,8 @@ CONTAINS
NULLIFY (qs_env%cp_ddapc_env)
NULLIFY (qs_env%cp_ddapc_ewald)
NULLIFY (qs_env%outer_scf_history)
NULLIFY (qs_env%gradient_history)
NULLIFY (qs_env%variable_history)
NULLIFY (qs_env%x_data)
NULLIFY (qs_env%et_coupling)
NULLIFY (qs_env%dftb_potential)
@ -965,6 +975,8 @@ CONTAINS
!> \param kpoints ...
!> \param WannierCentres ...
!> \param almo_scf_env ...
!> \param gradient_history ...
!> \param variable_history ...
!> \date 23.01.2002
!> \author MK
!> \version 1.0
@ -982,7 +994,7 @@ CONTAINS
outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, &
se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, ls_scf_env, &
do_transport, transport_env, lri_env, lri_density, dispersion_env, mp2_env, kg_env, &
kpoints, WannierCentres, almo_scf_env)
kpoints, WannierCentres, almo_scf_env, gradient_history, variable_history)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(cell_type), OPTIONAL, POINTER :: super_cell
@ -1048,6 +1060,7 @@ CONTAINS
TYPE(wannier_centres_type), DIMENSION(:), &
OPTIONAL, POINTER :: WannierCentres
TYPE(almo_scf_env_type), OPTIONAL, POINTER :: almo_scf_env
REAL(KIND=dp), DIMENSION(:, :), OPTIONAL, POINTER :: gradient_history, variable_history
CHARACTER(len=*), PARAMETER :: routineN = 'set_qs_env', routineP = moduleN//':'//routineN
@ -1074,6 +1087,8 @@ CONTAINS
! compiler (present in at least 10.0.025). A testcase has been submitted to intel.
IF (PRESENT(oce)) qs_env%oce => oce
IF (PRESENT(outer_scf_history)) qs_env%outer_scf_history => outer_scf_history
IF (PRESENT(gradient_history)) qs_env%gradient_history => gradient_history
IF (PRESENT(variable_history)) qs_env%variable_history => variable_history
IF (PRESENT(outer_scf_ihistory)) qs_env%outer_scf_ihistory = outer_scf_ihistory
IF (PRESENT(requires_mo_derivs)) qs_env%requires_mo_derivs = requires_mo_derivs
IF (PRESENT(has_unit_metric)) qs_env%has_unit_metric = has_unit_metric
@ -1310,6 +1325,10 @@ CONTAINS
DEALLOCATE (qs_env%outer_scf_history)
qs_env%outer_scf_ihistory = 0
ENDIF
IF (ASSOCIATED(qs_env%gradient_history)) &
DEALLOCATE (qs_env%gradient_history)
IF (ASSOCIATED(qs_env%variable_history)) &
DEALLOCATE (qs_env%variable_history)
IF (ASSOCIATED(qs_env%oce)) CALL deallocate_oce_set(qs_env%oce)
IF (ASSOCIATED(qs_env%local_rho_set)) THEN
CALL local_rho_set_release(qs_env%local_rho_set)

54
src/qs_ks_methods.F
View file

@ -686,34 +686,36 @@ CONTAINS
END IF
IF (calculate_forces .AND. dft_control%qs_control%becke_restraint) THEN
CALL pw_pool_give_back_pw(auxbas_pw_pool, becke%becke_pot%pw)
IF (dft_control%qs_control%becke_control%atomic_charges) THEN
DO iatom = 1, dft_control%qs_control%becke_control%natoms
CALL pw_pool_give_back_pw(auxbas_pw_pool, &
dft_control%qs_control%becke_control%charge(iatom)%pw)
END DO
END IF
IF (dft_control%qs_control%becke_control%cavity_confine) THEN
IF (.NOT. ASSOCIATED(becke%cavity_mat)) THEN
CALL pw_pool_give_back_pw(auxbas_pw_pool, becke%cavity%pw)
ELSE
DEALLOCATE (becke%cavity_mat)
IF (.NOT. becke%transfer_pot) THEN
CALL pw_pool_give_back_pw(auxbas_pw_pool, becke%becke_pot%pw)
IF (dft_control%qs_control%becke_control%atomic_charges) THEN
DO iatom = 1, dft_control%qs_control%becke_control%natoms
CALL pw_pool_give_back_pw(auxbas_pw_pool, &
dft_control%qs_control%becke_control%charge(iatom)%pw)
END DO
END IF
IF (dft_control%qs_control%becke_control%cavity_confine) THEN
IF (.NOT. ASSOCIATED(becke%cavity_mat)) THEN
CALL pw_pool_give_back_pw(auxbas_pw_pool, becke%cavity%pw)
ELSE
DEALLOCATE (becke%cavity_mat)
END IF
END IF
IF (ASSOCIATED(dft_control%qs_control%becke_control%charges_fragment)) &
DEALLOCATE (dft_control%qs_control%becke_control%charges_fragment)
IF (dft_control%qs_control%becke_control%constraint_type == cdft_combined_constraint) THEN
SELECT CASE (dft_control%qs_control%becke_control%combined_type)
CASE (cdft_combined_all)
! Do nothing
CASE (cdft_combined_acceptor, cdft_combined_donor)
CALL pw_pool_give_back_pw(auxbas_pw_pool, &
dft_control%qs_control%becke_control%combined_weight%pw)
END SELECT
END IF
dft_control%qs_control%becke_control%save_pot = .FALSE.
dft_control%qs_control%becke_control%need_pot = .TRUE.
dft_control%qs_control%becke_control%external_control = .FALSE.
END IF
IF (ASSOCIATED(dft_control%qs_control%becke_control%charges_fragment)) &
DEALLOCATE (dft_control%qs_control%becke_control%charges_fragment)
IF (dft_control%qs_control%becke_control%constraint_type == cdft_combined_constraint) THEN
SELECT CASE (dft_control%qs_control%becke_control%combined_type)
CASE (cdft_combined_all)
! Do nothing
CASE (cdft_combined_acceptor, cdft_combined_donor)
CALL pw_pool_give_back_pw(auxbas_pw_pool, &
dft_control%qs_control%becke_control%combined_weight%pw)
END SELECT
END IF
dft_control%qs_control%becke_control%save_pot = .FALSE.
dft_control%qs_control%becke_control%need_pot = .TRUE.
dft_control%qs_control%becke_control%external_control = .FALSE.
END IF
IF (calculate_forces .AND. dft_control%qs_control%cdft) THEN

174
src/qs_outer_scf.F
View file

@ -15,16 +15,19 @@ MODULE qs_outer_scf
ddapc_restraint_type,&
dft_control_type,&
s2_restraint_type
USE cp_log_handling, ONLY: cp_to_string
USE input_constants, ONLY: &
broyden_type_1, broyden_type_1_explicit, broyden_type_1_explicit_ls, broyden_type_1_ls, &
broyden_type_2, broyden_type_2_explicit, broyden_type_2_explicit_ls, broyden_type_2_ls, &
cdft2ot, cdft_combined_constraint, cdft_density_constraint, cdft_magnetization_constraint, &
do_ddapc_constraint, do_s2_constraint, ot2cdft, outer_scf_basis_center_opt, &
outer_scf_becke_constraint, outer_scf_cdft_constraint, outer_scf_ddapc_constraint, &
outer_scf_none, outer_scf_optimizer_bisect, outer_scf_optimizer_broyden, &
outer_scf_optimizer_diis, outer_scf_optimizer_newton, outer_scf_optimizer_none, &
outer_scf_optimizer_sd, outer_scf_s2_constraint
outer_scf_optimizer_diis, outer_scf_optimizer_newton, outer_scf_optimizer_newton_ls, &
outer_scf_optimizer_none, outer_scf_optimizer_sd, outer_scf_optimizer_secant, &
outer_scf_s2_constraint
USE kinds, ONLY: dp
USE mathlib, ONLY: diamat_all,&
invert_matrix
USE mathlib, ONLY: diamat_all
USE qs_basis_gradient, ONLY: qs_basis_center_gradient,&
qs_update_basis_center_pos,&
return_basis_center_gradient_norm
@ -203,9 +206,9 @@ CONTAINS
INTEGER :: handle, i, ibuf, ihigh, ihistory, ilow, &
j, jbuf, nb, nvar, optimizer_type
REAL(KIND=dp) :: interval, inv_error, scale, tmp
REAL(KIND=dp) :: interval, scale, tmp
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: ev
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: a, b, f, jacobian, x
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: a, b, f, x
REAL(KIND=dp), DIMENSION(:, :), POINTER :: inv_jacobian
CALL timeset(routineN, handle)
@ -301,66 +304,80 @@ CONTAINS
ENDDO
DEALLOCATE (a, b, ev)
ENDIF
CASE (outer_scf_optimizer_broyden)
CASE (outer_scf_optimizer_secant)
CPASSERT(SIZE(scf_env%outer_scf%gradient, 2) >= 3)
CPASSERT(SIZE(scf_env%outer_scf%gradient, 1) == 1)
nvar = SIZE(scf_env%outer_scf%gradient, 1)
IF (ihistory < 2) THEN
! Need two history values to use Broyden, switch to sd
! Need two history values to use secant, switch to sd
optimizer_type = outer_scf_optimizer_sd
CYCLE
END IF
IF (nvar == 1) THEN
! 1D case secant method
scf_env%outer_scf%variables(1, ihistory+1) = scf_env%outer_scf%variables(1, ihistory)- &
(scf_env%outer_scf%variables(1, ihistory)- &
scf_env%outer_scf%variables(1, ihistory-1))/ &
(scf_env%outer_scf%gradient(1, ihistory)- &
scf_env%outer_scf%gradient(1, ihistory-1))* &
scf_env%outer_scf%gradient(1, ihistory)
ELSE
! secant update
scf_env%outer_scf%variables(1, ihistory+1) = scf_env%outer_scf%variables(1, ihistory)- &
(scf_env%outer_scf%variables(1, ihistory)- &
scf_env%outer_scf%variables(1, ihistory-1))/ &
(scf_env%outer_scf%gradient(1, ihistory)- &
scf_env%outer_scf%gradient(1, ihistory-1))* &
scf_env%outer_scf%gradient(1, ihistory)
CASE (outer_scf_optimizer_broyden)
IF (.NOT. ASSOCIATED(scf_env%outer_scf%inv_jacobian)) THEN
! Inverse Jacobian not yet built, switch to sd
optimizer_type = outer_scf_optimizer_sd
CYCLE
END IF
inv_jacobian => scf_env%outer_scf%inv_jacobian
IF (ihistory < 2) THEN
! Cannot perform a Broyden update without enough SCF history on this energy evaluation
scf_control%outer_scf%broyden_update = .FALSE.
END IF
IF (scf_control%outer_scf%broyden_update) THEN
! Perform a Broyden update of the inverse Jacobian J^(-1)
IF (SIZE(scf_env%outer_scf%gradient, 2) .LT. 3) &
CALL cp_abort(__LOCATION__, &
"Keyword EXTRAPOLATION_ORDER in section OUTER_SCF "// &
"must be greater than or equal to 3 for Broyden optimizers.")
nvar = SIZE(scf_env%outer_scf%gradient, 1)
ALLOCATE (f(nvar, 1), x(nvar, 1))
DO i = 1, nvar
f(i, 1) = scf_env%outer_scf%gradient(i, ihistory)-scf_env%outer_scf%gradient(i, ihistory-1)
x(i, 1) = scf_env%outer_scf%variables(i, ihistory)-scf_env%outer_scf%variables(i, ihistory-1)
END DO
! Use finite difference Jacobian as initial guess
! TODO: Add possibility to restart Jacobian
IF (ihistory == 2) THEN
IF (.NOT. ASSOCIATED(scf_env%outer_scf%inv_jacobian)) &
ALLOCATE (scf_env%outer_scf%inv_jacobian(nvar, nvar))
inv_jacobian => scf_env%outer_scf%inv_jacobian
ALLOCATE (jacobian(nvar, nvar))
jacobian = 0.0_dp
DO i = 1, nvar
jacobian(i, i) = f(i, 1)/x(i, 1)
END DO
CALL invert_matrix(jacobian, inv_jacobian, inv_error)
SELECT CASE (scf_control%outer_scf%broyden_type)
CASE (broyden_type_1, broyden_type_1_explicit, broyden_type_1_ls, broyden_type_1_explicit_ls)
! Broyden's 1st method
! Denote: dx_n = \delta x_n; df_n = \delta f_n
! J_(n+1)^(-1) = J_n^(-1) + (dx_n - J_n^(-1)*df_n)*(dx_n^T * J_n^(-1))/(dx_n^T * J_n^(-1) * df_n)
scale = SUM(MATMUL(TRANSPOSE(x), MATMUL(inv_jacobian, f)))
DEALLOCATE (jacobian)
ELSE
SELECT CASE (scf_control%outer_scf%jacobian_type)
CASE DEFAULT
! Broyden's 1st method
! Update inverse Jacobian: dx_n = \delta x_n; df_n = \delta f_n
! J_(n+1)^(-1) = J_n^(-1) + (dx_n - J_n^(-1)*df_n)*(dx_n^T * J_n^(-1))/(dx_n^T * J_n^(-1) * df_n)
inv_jacobian => scf_env%outer_scf%inv_jacobian
scale = SUM(MATMUL(TRANSPOSE(x), MATMUL(inv_jacobian, f)))
scale = 1.0_dp/scale
IF (scale < 1.0E-12_dp) scale = 1.0E-12_dp
inv_jacobian = inv_jacobian+scale*MATMUL((x-MATMUL(inv_jacobian, f)), &
MATMUL(TRANSPOSE(x), inv_jacobian))
END SELECT
END IF
! Broyden update: x_(n+1) = x_n - J^(-1)*f(x_n)
scf_env%outer_scf%variables(:, ihistory+1) = scf_env%outer_scf%variables(:, ihistory)- &
MATMUL(inv_jacobian, scf_env%outer_scf%gradient(:, ihistory))
scale = 1.0_dp/scale
IF (scale < 1.0E-12_dp) scale = 1.0E-12_dp
inv_jacobian = inv_jacobian+scale*MATMUL((x-MATMUL(inv_jacobian, f)), &
MATMUL(TRANSPOSE(x), inv_jacobian))
CASE (broyden_type_2, broyden_type_2_explicit, broyden_type_2_ls, broyden_type_2_explicit_ls)
! Broyden's 2nd method
! J_(n+1)^(-1) = J_n^(-1) + (dx_n - J_n^(-1)*df_n)*(df_n^T)/(||df_n||^2)
scale = SUM(MATMUL(TRANSPOSE(f), f))
scale = 1.0_dp/scale
IF (scale < 1.0E-12_dp) scale = 1.0E-12_dp
inv_jacobian = inv_jacobian+scale*MATMUL((x-MATMUL(inv_jacobian, f)), TRANSPOSE(inv_jacobian))
CASE DEFAULT
CALL cp_abort(__LOCATION__, &
"Unknown Broyden type: "// &
cp_to_string(scf_control%outer_scf%broyden_type))
END SELECT
! Clean up
DEALLOCATE (f, x)
END IF
CASE (outer_scf_optimizer_newton)
! Update variables x_(n+1) = x_n - J^(-1)*f(x_n)
scf_env%outer_scf%variables(:, ihistory+1) = scf_env%outer_scf%variables(:, ihistory)- &
scf_control%outer_scf%newton_step* &
MATMUL(inv_jacobian, scf_env%outer_scf%gradient(:, ihistory))
scf_control%outer_scf%broyden_update = .TRUE.
CASE (outer_scf_optimizer_newton, outer_scf_optimizer_newton_ls)
CPASSERT(ASSOCIATED(scf_env%outer_scf%inv_jacobian))
inv_jacobian => scf_env%outer_scf%inv_jacobian
scf_env%outer_scf%variables(:, ihistory+1) = scf_env%outer_scf%variables(:, ihistory)- &
scf_control%outer_scf%newton_step* &
MATMUL(inv_jacobian, scf_env%outer_scf%gradient(:, ihistory))
CASE DEFAULT
CPABORT("")
@ -501,7 +518,6 @@ CONTAINS
CPABORT("")
END SELECT
CALL set_qs_env(qs_env, outer_scf_ihistory=outer_scf_ihistory)
! multilinear extrapolation
nvec = MIN(nhistory, outer_scf_ihistory)
alpha = nvec
@ -558,10 +574,16 @@ CONTAINS
CHARACTER(len=*), PARAMETER :: routineN = 'outer_loop_switch', &
routineP = moduleN//':'//routineN
INTEGER :: nvariables
SELECT CASE (dir)
! Constraint -> OT
CASE (cdft2ot)
! Switch data in scf_control
! Constraint -> OT
! Switch data in scf_control: first save values that might have changed
cdft_control%constraint_control%broyden_update = scf_control%outer_scf%broyden_update
cdft_control%constraint_control%newton_step = scf_control%outer_scf%newton_step
cdft_control%constraint_control%ijacobian = scf_control%outer_scf%ijacobian
! Now switch
scf_control%outer_scf%have_scf = cdft_control%ot_control%have_scf
scf_control%outer_scf%max_scf = cdft_control%ot_control%max_scf
scf_control%outer_scf%eps_scf = cdft_control%ot_control%eps_scf
@ -571,12 +593,26 @@ CONTAINS
scf_control%outer_scf%diis_buffer_length = cdft_control%ot_control%diis_buffer_length
scf_control%outer_scf%bisect_trust_count = cdft_control%ot_control%bisect_trust_count
scf_control%outer_scf%jacobian_type = cdft_control%ot_control%jacobian_type
scf_control%outer_scf%jacobian_step = cdft_control%ot_control%jacobian_step
scf_control%outer_scf%jacobian_freq = cdft_control%ot_control%jacobian_freq
scf_control%outer_scf%newton_step = cdft_control%ot_control%newton_step
scf_control%outer_scf%ijacobian = cdft_control%ot_control%ijacobian
scf_control%outer_scf%build_jacobian = cdft_control%ot_control%build_jacobian
scf_control%outer_scf%broyden_type = cdft_control%ot_control%broyden_type
scf_control%outer_scf%broyden_update = cdft_control%ot_control%broyden_update
! Switch data in scf_env: first save current values for constraint
cdft_control%constraint%iter_count = scf_env%outer_scf%iter_count
cdft_control%constraint%energy = scf_env%outer_scf%energy
cdft_control%constraint%variables = scf_env%outer_scf%variables
cdft_control%constraint%gradient = scf_env%outer_scf%gradient
cdft_control%constraint%count = scf_env%outer_scf%count
cdft_control%constraint%deallocate_jacobian = scf_env%outer_scf%deallocate_jacobian
IF (ASSOCIATED(scf_env%outer_scf%inv_jacobian)) THEN
nvariables = SIZE(scf_env%outer_scf%inv_jacobian, 1)
IF (.NOT. ASSOCIATED(cdft_control%constraint%inv_jacobian)) &
ALLOCATE (cdft_control%constraint%inv_jacobian(nvariables, nvariables))
cdft_control%constraint%inv_jacobian = scf_env%outer_scf%inv_jacobian
END IF
! Now switch
IF (ASSOCIATED(scf_env%outer_scf%energy)) &
DEALLOCATE (scf_env%outer_scf%energy)
@ -594,8 +630,12 @@ CONTAINS
DEALLOCATE (scf_env%outer_scf%count)
ALLOCATE (scf_env%outer_scf%count(scf_control%outer_scf%max_scf+1))
scf_env%outer_scf%count = 0
! OT -> constraint
! OT SCF does not need Jacobian
scf_env%outer_scf%deallocate_jacobian = .TRUE.
IF (ASSOCIATED(scf_env%outer_scf%inv_jacobian)) &
DEALLOCATE (scf_env%outer_scf%inv_jacobian)
CASE (ot2cdft)
! OT -> constraint
scf_control%outer_scf%have_scf = cdft_control%constraint_control%have_scf
scf_control%outer_scf%max_scf = cdft_control%constraint_control%max_scf
scf_control%outer_scf%eps_scf = cdft_control%constraint_control%eps_scf
@ -605,25 +645,38 @@ CONTAINS
scf_control%outer_scf%diis_buffer_length = cdft_control%constraint_control%diis_buffer_length
scf_control%outer_scf%bisect_trust_count = cdft_control%constraint_control%bisect_trust_count
scf_control%outer_scf%jacobian_type = cdft_control%constraint_control%jacobian_type
scf_control%outer_scf%jacobian_step = cdft_control%constraint_control%jacobian_step
scf_control%outer_scf%jacobian_freq = cdft_control%constraint_control%jacobian_freq
scf_control%outer_scf%ijacobian = cdft_control%constraint_control%ijacobian
scf_control%outer_scf%build_jacobian = cdft_control%constraint_control%build_jacobian
scf_control%outer_scf%broyden_type = cdft_control%constraint_control%broyden_type
scf_control%outer_scf%broyden_update = cdft_control%constraint_control%broyden_update
scf_control%outer_scf%newton_step = cdft_control%constraint_control%newton_step
nvariables = SIZE(cdft_control%constraint%variables, 1)
IF (ASSOCIATED(scf_env%outer_scf%energy)) &
DEALLOCATE (scf_env%outer_scf%energy)
ALLOCATE (scf_env%outer_scf%energy(scf_control%outer_scf%max_scf+1))
scf_env%outer_scf%energy = cdft_control%constraint%energy
IF (ASSOCIATED(scf_env%outer_scf%variables)) &
DEALLOCATE (scf_env%outer_scf%variables)
ALLOCATE (scf_env%outer_scf%variables(SIZE(cdft_control%constraint%variables, 1), &
scf_control%outer_scf%max_scf+1))
ALLOCATE (scf_env%outer_scf%variables(nvariables, scf_control%outer_scf%max_scf+1))
scf_env%outer_scf%variables = cdft_control%constraint%variables
IF (ASSOCIATED(scf_env%outer_scf%gradient)) &
DEALLOCATE (scf_env%outer_scf%gradient)
ALLOCATE (scf_env%outer_scf%gradient(SIZE(cdft_control%constraint%gradient, 1), &
scf_control%outer_scf%max_scf+1))
ALLOCATE (scf_env%outer_scf%gradient(nvariables, scf_control%outer_scf%max_scf+1))
scf_env%outer_scf%gradient = cdft_control%constraint%gradient
IF (ASSOCIATED(scf_env%outer_scf%count)) &
DEALLOCATE (scf_env%outer_scf%count)
ALLOCATE (scf_env%outer_scf%count(scf_control%outer_scf%max_scf+1))
scf_env%outer_scf%count = cdft_control%constraint%count
scf_env%outer_scf%iter_count = cdft_control%constraint%iter_count
scf_env%outer_scf%deallocate_jacobian = cdft_control%constraint%deallocate_jacobian
IF (ASSOCIATED(cdft_control%constraint%inv_jacobian)) THEN
IF (ASSOCIATED(scf_env%outer_scf%inv_jacobian)) &
DEALLOCATE (scf_env%outer_scf%inv_jacobian)
ALLOCATE (scf_env%outer_scf%inv_jacobian(nvariables, nvariables))
scf_env%outer_scf%inv_jacobian = cdft_control%constraint%inv_jacobian
END IF
CASE DEFAULT
CPABORT("")
END SELECT
@ -644,15 +697,20 @@ CONTAINS
routineP = moduleN//':'//routineN
INTEGER :: handle, outer_scf_ihistory
REAL(kind=dp), DIMENSION(:, :), POINTER :: outer_scf_history
REAL(kind=dp), DIMENSION(:, :), POINTER :: gradient_history, outer_scf_history, &
variable_history
CALL timeset(routineN, handle)
CALL get_qs_env(qs_env, outer_scf_history=outer_scf_history, &
outer_scf_ihistory=outer_scf_ihistory)
outer_scf_ihistory=outer_scf_ihistory, &
gradient_history=gradient_history, &
variable_history=variable_history)
CPASSERT(SIZE(outer_scf_history, 2) > 0)
outer_scf_ihistory = 0
outer_scf_history = 0.0_dp
gradient_history = 0.0_dp
variable_history = 0.0_dp
CALL set_qs_env(qs_env, outer_scf_ihistory=outer_scf_ihistory)
CALL timestop(handle)

677
src/qs_scf.F
View file

@ -46,14 +46,18 @@ MODULE qs_scf
USE atomic_kind_types, ONLY: atomic_kind_type
USE cp_control_types, ONLY: cdft_control_type,&
dft_control_type
USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
copy_fm_to_dbcsr
USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm
USE cp_files, ONLY: open_file
USE cp_fm_types, ONLY: cp_fm_create,&
cp_fm_release,&
cp_fm_to_fm,&
cp_fm_type
USE cp_log_handling, ONLY: cp_get_default_logger,&
USE cp_log_handling, ONLY: cp_add_default_logger,&
cp_get_default_logger,&
cp_logger_create,&
cp_logger_release,&
cp_logger_type,&
cp_rm_default_logger,&
cp_to_string
USE cp_output_handling, ONLY: cp_add_iter_level,&
cp_iterate,&
@ -74,15 +78,18 @@ MODULE qs_scf
dbcsr_set,&
dbcsr_type
USE input_constants, ONLY: &
broyden_type_1, broyden_type_1_explicit, broyden_type_1_explicit_ls, broyden_type_1_ls, &
broyden_type_2, broyden_type_2_explicit, broyden_type_2_explicit_ls, broyden_type_2_ls, &
cdft2ot, cdft_combined_acceptor, cdft_combined_all, cdft_combined_constraint, &
cdft_combined_donor, history_guess, jacobian_fd1, jacobian_fd1_backward, &
jacobian_fd1_central, jacobian_fd2, jacobian_fd2_backward, ot2cdft, ot_precond_full_all, &
ot_precond_full_single, ot_precond_full_single_inverse, ot_precond_none, &
ot_precond_s_inverse, outer_scf_hirshfeld_constraint, outer_scf_optimizer_broyden, &
outer_scf_optimizer_newton
outer_scf_optimizer_newton, outer_scf_optimizer_newton_ls
USE input_section_types, ONLY: section_vals_get_subs_vals,&
section_vals_type
USE kinds, ONLY: default_string_length,&
USE kinds, ONLY: default_path_length,&
default_string_length,&
dp
USE kpoint_io, ONLY: write_kpoints_restart
USE kpoint_types, ONLY: kpoint_type
@ -162,6 +169,7 @@ MODULE qs_scf
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_scf'
LOGICAL, PRIVATE :: reuse_precond = .FALSE.
LOGICAL, PRIVATE :: used_history = .FALSE.
PUBLIC :: scf, scf_env_cleanup, scf_env_do_scf, cdft_scf
@ -180,9 +188,10 @@ CONTAINS
CHARACTER(len=*), PARAMETER :: routineN = 'scf', routineP = moduleN//':'//routineN
INTEGER :: max_scf_tmp
INTEGER :: ihistory, max_scf_tmp
LOGICAL :: converged, outer_scf_loop, should_stop
LOGICAL, SAVE :: first_step_flag = .TRUE.
REAL(KIND=dp), DIMENSION(:, :), POINTER :: gradient_history, variable_history
TYPE(cp_logger_type), POINTER :: logger
TYPE(dft_control_type), POINTER :: dft_control
TYPE(qs_scf_env_type), POINTER :: scf_env
@ -227,6 +236,25 @@ CONTAINS
! If SCF has not converged, then we should not start MP2
IF (ASSOCIATED(qs_env%mp2_env)) qs_env%mp2_env%hf_fail = .NOT. converged
! Add the converged outer_scf SCF gradient(s)/variable(s) to history
IF (scf_control%outer_scf%have_scf) THEN
ihistory = scf_env%outer_scf%iter_count
CALL get_qs_env(qs_env, gradient_history=gradient_history, &
variable_history=variable_history)
! We only store the latest two values
gradient_history(:, 1) = gradient_history(:, 2)
gradient_history(:, 2) = scf_env%outer_scf%gradient(:, ihistory)
variable_history(:, 1) = variable_history(:, 2)
variable_history(:, 2) = scf_env%outer_scf%variables(:, ihistory)
! Reset flag
IF (used_history) used_history = .FALSE.
! Update a counter and check if the Jacobian should be deallocated
IF (ASSOCIATED(scf_env%outer_scf%inv_jacobian)) THEN
scf_control%outer_scf%ijacobian = scf_control%outer_scf%ijacobian+1
IF (scf_control%outer_scf%ijacobian .GE. scf_control%outer_scf%jacobian_freq(2)) &
scf_env%outer_scf%deallocate_jacobian = .TRUE.
END IF
END IF
! *** add the converged wavefunction to the wavefunction history
IF ((ASSOCIATED(qs_env%wf_history)) .AND. &
((scf_control%density_guess .NE. history_guess) .OR. &
@ -494,6 +522,9 @@ CONTAINS
IF (exit_outer_loop) EXIT scf_outer_loop
! Check if the inverse Jacobian needs to be calculated
CALL qs_calculate_inverse_jacobian(qs_env)
!
CALL outer_loop_optimize(scf_env, scf_control)
CALL outer_loop_update_qs_env(qs_env, scf_env)
CALL qs_ks_did_change(ks_env, potential_changed=.TRUE.)
@ -879,7 +910,8 @@ CONTAINS
IF (ASSOCIATED(scf_env%outer_scf%energy)) THEN
DEALLOCATE (scf_env%outer_scf%energy)
ENDIF
IF (ASSOCIATED(scf_env%outer_scf%inv_jacobian)) THEN
IF (ASSOCIATED(scf_env%outer_scf%inv_jacobian) .AND. &
scf_env%outer_scf%deallocate_jacobian) THEN
DEALLOCATE (scf_env%outer_scf%inv_jacobian)
ENDIF
@ -995,14 +1027,18 @@ CONTAINS
energy, cdft_control%total_steps, &
should_stop, cdft_loop_converged, cdft_loop=.TRUE.)
IF (exit_cdft_loop) EXIT cdft_outer_loop
! Check if inverse Jacobian needs to be calculated
CALL outer_loop_calculate_inverse_jacobian(qs_env)
! Check if the inverse Jacobian needs to be calculated
CALL qs_calculate_inverse_jacobian(qs_env)
! Check if a line search should be performed to find an optimal step size for the optimizer
CALL qs_cdft_line_search(qs_env)
! Optimize constraint
CALL outer_loop_optimize(scf_env, scf_control)
CALL outer_loop_update_qs_env(qs_env, scf_env)
CALL qs_ks_did_change(ks_env, potential_changed=.TRUE.)
END DO cdft_outer_loop
cdft_control%ienergy = cdft_control%ienergy+1
IF (dft_control%qs_control%becke_restraint) THEN
! Store needed arrays for ET coupling calculation
IF (cdft_control%do_et) THEN
@ -1040,8 +1076,9 @@ CONTAINS
END DO
END IF
END IF
! Forces are not needed => time to deallocate Becke
IF (.NOT. dft_control%qs_control%becke_control%save_pot) THEN
! Forces are not needed/serial mode if mixed CDFT calculation is not active => time to deallocate Becke
IF (.NOT. (dft_control%qs_control%becke_control%save_pot .OR. &
dft_control%qs_control%becke_control%transfer_pot)) THEN
CALL get_qs_env(qs_env, pw_env=pw_env)
CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
CALL pw_pool_give_back_pw(auxbas_pw_pool, &
@ -1112,6 +1149,9 @@ CONTAINS
DEALLOCATE (cdft_control%constraint%gradient)
IF (ASSOCIATED(cdft_control%constraint%energy)) &
DEALLOCATE (cdft_control%constraint%energy)
IF (ASSOCIATED(cdft_control%constraint%inv_jacobian) .AND. &
cdft_control%constraint%deallocate_jacobian) &
DEALLOCATE (cdft_control%constraint%inv_jacobian)
END SUBROUTINE cdft_control_cleanup
@ -1121,60 +1161,434 @@ CONTAINS
!> \par History
!> 01.2017 created [Nico Holmberg]
! **************************************************************************************************
SUBROUTINE outer_loop_calculate_inverse_jacobian(qs_env)
SUBROUTINE qs_calculate_inverse_jacobian(qs_env)
TYPE(qs_environment_type), POINTER :: qs_env
CHARACTER(LEN=*), PARAMETER :: routineN = 'outer_loop_calculate_inverse_jacobian', &
CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_calculate_inverse_jacobian', &
routineP = moduleN//':'//routineN
INTEGER :: handle, i, ispin, iter_count, iwork, j, &
max_scf, nspins, nvar, nwork, &
CHARACTER(len=default_path_length) :: project_name
INTEGER :: handle, i, ihistory, ispin, iter_count, &
iwork, j, lp, max_scf, nspins, nvar, &
nwork, outer_scf_ihistory, &
output_unit, pwork
LOGICAL :: converged, should_stop
LOGICAL :: converged, explicit_jacobian, &
should_build, should_stop, &
use_md_history
REAL(KIND=dp) :: dh, inv_error
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: coeff, step_multiplier
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: jacobian
REAL(KIND=dp), DIMENSION(:), POINTER :: energy
REAL(KIND=dp), DIMENSION(:, :), POINTER :: gradient, inv_jacobian
REAL(KIND=dp), DIMENSION(:, :), POINTER :: gradient, gradient_history, &
inv_jacobian, variable_history
TYPE(cdft_control_type), POINTER :: cdft_control
TYPE(cp_logger_type), POINTER :: logger
TYPE(cp_logger_type), POINTER :: logger, tmp_logger
TYPE(cp_para_env_type), POINTER :: para_env
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: p_rmpv
TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao_kp
TYPE(dft_control_type), POINTER :: dft_control
TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mos, mos_stashed
TYPE(qs_energy_type), POINTER :: energy_qs
TYPE(qs_ks_env_type), POINTER :: ks_env
TYPE(qs_rho_type), POINTER :: rho
TYPE(qs_scf_env_type), POINTER :: scf_env
TYPE(scf_control_type), POINTER :: scf_control
TYPE(section_vals_type), POINTER :: input, scf_section
CALL timeset(routineN, handle)
NULLIFY (energy, gradient, p_rmpv, rho_ao_kp, mos, rho, &
mos_stashed, ks_env, scf_env, scf_control, dft_control, &
cdft_control, input, scf_section, inv_jacobian)
cdft_control, inv_jacobian, para_env, &
tmp_logger, gradient_history, variable_history, energy_qs)
logger => cp_get_default_logger()
CPASSERT(ASSOCIATED(qs_env))
CALL get_qs_env(qs_env, scf_env=scf_env, ks_env=ks_env, &
scf_control=scf_control, mos=mos, rho=rho, &
dft_control=dft_control, input=input)
dft_control=dft_control, &
para_env=para_env, energy=energy_qs)
explicit_jacobian = .FALSE.
should_build = .FALSE.
use_md_history = .FALSE.
iter_count = scf_env%outer_scf%iter_count
SELECT CASE (scf_control%outer_scf%optimizer)
CASE DEFAULT
scf_control%outer_scf%build_jacobian = .FALSE.
CASE (outer_scf_optimizer_newton)
CASE (outer_scf_optimizer_newton, outer_scf_optimizer_newton_ls)
scf_control%outer_scf%build_jacobian = .TRUE.
explicit_jacobian = .TRUE.
CASE (outer_scf_optimizer_broyden)
! TODO: Add Broyden variants where the explicit Jacobian is built
scf_control%outer_scf%build_jacobian = .FALSE.
SELECT CASE (scf_control%outer_scf%broyden_type)
CASE (broyden_type_1, broyden_type_2, broyden_type_1_ls, broyden_type_2_ls)
scf_control%outer_scf%build_jacobian = .TRUE.
explicit_jacobian = .FALSE.
CASE (broyden_type_1_explicit, broyden_type_2_explicit, broyden_type_1_explicit_ls, broyden_type_2_explicit_ls)
scf_control%outer_scf%build_jacobian = .TRUE.
explicit_jacobian = .TRUE.
END SELECT
END SELECT
IF (scf_control%outer_scf%build_jacobian) THEN
! Get output_unit
scf_section => section_vals_get_subs_vals(input, "DFT%SCF")
output_unit = cp_print_key_unit_nr(logger, scf_section, "PRINT%PROGRAM_RUN_INFO", &
extension=".scfLog")
! Check if an old Jacobian can be reused avoiding a rebuild
IF (ASSOCIATED(scf_env%outer_scf%inv_jacobian)) THEN
! Rebuild if number of previous energy evaluations exceeds limit
IF (scf_control%outer_scf%ijacobian .GE. scf_control%outer_scf%jacobian_freq(2) &
.AND. .NOT. used_history) THEN
should_build = .TRUE.
! Zero the corresponding counter
scf_control%outer_scf%ijacobian = 0
! Rebuild if number of previous SCF iterations exceeds limit (on this energy eval)
ELSE IF (MODULO(iter_count, scf_control%outer_scf%jacobian_freq(1)) == 0) THEN
should_build = .TRUE.
END IF
IF (should_build) DEALLOCATE (scf_env%outer_scf%inv_jacobian)
ELSE
should_build = .TRUE.
! Zero the counter
scf_control%outer_scf%ijacobian = 0
END IF
END IF
IF (should_build) THEN
scf_env%outer_scf%deallocate_jacobian = .FALSE.
! Actually need to (re)build the Jacobian
IF (explicit_jacobian) THEN
! Build Jacobian with finite differences
cdft_control => dft_control%qs_control%cdft_control
IF (.NOT. ASSOCIATED(cdft_control)) &
CALL cp_abort(__LOCATION__, &
"Optimizers that need the explicit Jacobian can"// &
" only be used together with a valid CDFT constraint")
! Redirect output from Jacobian calculation to a new file by creating a temporary logger
project_name = logger%iter_info%project_name
IF (para_env%mepos == para_env%source) THEN
lp = LEN_TRIM(project_name)
project_name(lp+1:LEN(project_name)) = "-JacobianInfo.out"
CALL open_file(file_name=project_name, file_status="UNKNOWN", &
file_action="WRITE", file_position="APPEND", &
unit_number=output_unit)
ELSE
output_unit = -1
END IF
CALL cp_logger_create(tmp_logger, &
para_env=para_env, &
default_global_unit_nr=output_unit, &
close_global_unit_on_dealloc=.FALSE.)
! Save last converged state so we can roll back to it (mo_coeff and some outer_loop variables)
nspins = dft_control%nspins
ALLOCATE (mos_stashed(nspins))
DO ispin = 1, nspins
CALL duplicate_mo_set(mos_stashed(ispin)%mo_set, mos(ispin)%mo_set)
END DO
CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
p_rmpv => rho_ao_kp(:, 1)
! Allocate work
nvar = SIZE(scf_env%outer_scf%variables, 1)
max_scf = scf_control%outer_scf%max_scf+1
ALLOCATE (gradient(nvar, max_scf))
gradient = scf_env%outer_scf%gradient
ALLOCATE (energy(max_scf))
energy = scf_env%outer_scf%energy
ALLOCATE (jacobian(nvar, nvar))
jacobian = 0.0_dp
! Setup finite difference scheme
SELECT CASE (scf_control%outer_scf%jacobian_type)
CASE DEFAULT
CALL cp_abort(__LOCATION__, &
"Unknown Jacobian type: "// &
cp_to_string(scf_control%outer_scf%jacobian_type))
CASE (jacobian_fd1)
! f'(x0) = [ f(x0+h) - f(x0) ] / h
nwork = 0
pwork = 1
ALLOCATE (coeff(nwork:pwork), step_multiplier(nwork:pwork))
coeff(nwork) = -1.0_dp
coeff(pwork) = 1.0_dp
step_multiplier = 1.0_dp
dh = scf_control%outer_scf%jacobian_step
CASE (jacobian_fd1_backward)
! f'(x0) = [ f(x0) - f(x0-h) ] / h
nwork = -1
pwork = 0
ALLOCATE (coeff(nwork:pwork), step_multiplier(nwork:pwork))
coeff(nwork) = -1.0_dp
coeff(pwork) = 1.0_dp
step_multiplier = -1.0_dp
dh = scf_control%outer_scf%jacobian_step
CASE (jacobian_fd2)
! f'(x0) = [ -f(x0+2h) + 4f(x0+h) - 3f(x0) ] / 2h
nwork = 0
pwork = 2
ALLOCATE (coeff(nwork:pwork), step_multiplier(nwork:pwork))
coeff(0) = -3.0_dp
coeff(1) = 4.0_dp
coeff(2) = -1.0_dp
step_multiplier(0) = 0.0_dp
step_multiplier(1) = 1.0_dp
step_multiplier(2) = 2.0_dp
dh = 2.0_dp*scf_control%outer_scf%jacobian_step
CASE (jacobian_fd2_backward)
! f'(x0) = [ 3f(x0) - 4f(x0-h) + f(x0-2h) ] / 2h
nwork = -2
pwork = 0
ALLOCATE (coeff(nwork:pwork), step_multiplier(nwork:pwork))
coeff(0) = 3.0_dp
coeff(-1) = -4.0_dp
coeff(-2) = 1.0_dp
step_multiplier(0) = 0.0_dp
step_multiplier(-1) = -1.0_dp
step_multiplier(-2) = -2.0_dp
dh = 2.0_dp*scf_control%outer_scf%jacobian_step
CASE (jacobian_fd1_central)
! f'(x0) = [ f(x0+h) - f(x0-h) ] / 2h
nwork = -1
pwork = 1
ALLOCATE (coeff(nwork:pwork), step_multiplier(nwork:pwork))
coeff(0) = 0.0_dp
coeff(nwork) = -1.0_dp
coeff(pwork) = 1.0_dp
step_multiplier(0) = 0.0_dp
step_multiplier(nwork) = -1.0_dp
step_multiplier(pwork) = 1.0_dp
dh = 2.0_dp*scf_control%outer_scf%jacobian_step
END SELECT
DO i = 1, nvar
jacobian(i, :) = coeff(0)*scf_env%outer_scf%gradient(i, iter_count)
END DO
! Print some info
IF (output_unit > 0) THEN
WRITE (output_unit, FMT="(/,A)") &
" ================================== JACOBIAN CALCULATION ================================="
WRITE (output_unit, FMT="(A)") &
" Evaluating inverse Jacobian using finite differences"
WRITE (output_unit, '(A,I10,A,I10)') &
" Energy evaluation: ", cdft_control%ienergy, ", CDFT SCF iteration: ", iter_count
SELECT CASE (scf_control%outer_scf%jacobian_type)
CASE (jacobian_fd1)
WRITE (output_unit, '(A)') " Type : First order forward difference"
CASE (jacobian_fd1_backward)
WRITE (output_unit, '(A)') " Type : First order backward difference"
CASE (jacobian_fd2)
WRITE (output_unit, '(A)') " Type : Second order forward difference"
CASE (jacobian_fd2_backward)
WRITE (output_unit, '(A)') " Type : Second order backward difference"
CASE (jacobian_fd1_central)
WRITE (output_unit, '(A)') " Type : First order central difference"
CASE DEFAULT
CALL cp_abort(__LOCATION__, "Unknown Jacobian type: "// &
cp_to_string(scf_control%outer_scf%jacobian_type))
END SELECT
WRITE (output_unit, '(A,ES12.4)') " Step size : ", scf_control%outer_scf%jacobian_step
END IF
! Calculate the Jacobian by perturbing each Lagrangian and recalculating the energy self-consistently
CALL cp_add_default_logger(tmp_logger)
DO i = 1, nvar
DO iwork = nwork, pwork
IF (iwork == 0) CYCLE
scf_env%outer_scf%variables(i, iter_count+1) = scf_env%outer_scf%variables(i, iter_count)+ &
step_multiplier(iwork)* &
scf_control%outer_scf%jacobian_step
CALL outer_loop_update_qs_env(qs_env, scf_env)
CALL qs_ks_did_change(ks_env, potential_changed=.TRUE.)
CALL outer_loop_switch(scf_env, scf_control, cdft_control, cdft2ot)
CALL scf_env_do_scf(scf_env=scf_env, scf_control=scf_control, qs_env=qs_env, &
converged=converged, should_stop=should_stop)
CALL outer_loop_switch(scf_env, scf_control, cdft_control, ot2cdft)
! Update (iter_count + 1) element of gradient and print constraint info
scf_env%outer_scf%iter_count = scf_env%outer_scf%iter_count+1
CALL outer_loop_gradient(qs_env, scf_env)
CALL qs_scf_cdft_info(output_unit, scf_control, scf_env, cdft_control, &
energy_qs, cdft_control%total_steps, &
should_stop=.FALSE., outer_loop_converged=.FALSE., cdft_loop=.FALSE.)
scf_env%outer_scf%iter_count = scf_env%outer_scf%iter_count-1
! Update Jacobian
DO j = 1, nvar
jacobian(j, i) = jacobian(j, i)+coeff(iwork)*scf_env%outer_scf%gradient(j, iter_count+1)
END DO
! Reset everything to last converged state
scf_env%outer_scf%variables(i, iter_count+1) = 0.0_dp
scf_env%outer_scf%gradient = gradient
scf_env%outer_scf%energy = energy
DO ispin = 1, nspins
CALL deallocate_mo_set(mos(ispin)%mo_set)
CALL duplicate_mo_set(mos(ispin)%mo_set, mos_stashed(ispin)%mo_set)
CALL calculate_density_matrix(mos(ispin)%mo_set, &
p_rmpv(ispin)%matrix)
END DO
CALL qs_rho_update_rho(rho, qs_env=qs_env)
CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.)
END DO
END DO
CALL cp_rm_default_logger()
CALL cp_logger_release(tmp_logger)
! Finalize and invert Jacobian
jacobian = jacobian/dh
IF (.NOT. ASSOCIATED(scf_env%outer_scf%inv_jacobian)) &
ALLOCATE (scf_env%outer_scf%inv_jacobian(nvar, nvar))
inv_jacobian => scf_env%outer_scf%inv_jacobian
CALL invert_matrix(jacobian, inv_jacobian, inv_error)
scf_control%outer_scf%broyden_update = .FALSE.
! Release temporary storage
DO ispin = 1, nspins
CALL deallocate_mo_set(mos_stashed(ispin)%mo_set)
END DO
DEALLOCATE (mos_stashed, jacobian, gradient, energy, coeff, step_multiplier)
IF (output_unit > 0) THEN
WRITE (output_unit, FMT="(/,A)") &
" ================================== JACOBIAN CALCULATED =================================="
END IF
ELSE
! Build a strictly diagonal Jacobian from history and invert it
CALL get_qs_env(qs_env, outer_scf_ihistory=outer_scf_ihistory)
ihistory = scf_env%outer_scf%iter_count
IF (outer_scf_ihistory .GE. 3 .AND. .NOT. used_history) THEN
! First, lets try using the history from previous energy evaluations
CALL get_qs_env(qs_env, gradient_history=gradient_history, &
variable_history=variable_history)
nvar = SIZE(scf_env%outer_scf%variables, 1)
use_md_history = .TRUE.
! Check that none of the history values are identical in which case we should try something different
DO i = 1, nvar
IF (ABS(variable_history(i, 2)-variable_history(i, 1)) .LT. 1.0E-12_dp) &
use_md_history = .FALSE.
END DO
IF (use_md_history) THEN
ALLOCATE (jacobian(nvar, nvar))
DO i = 1, nvar
jacobian(i, i) = (gradient_history(i, 2)-gradient_history(i, 1))/ &
(variable_history(i, 2)-variable_history(i, 1))
END DO
IF (.NOT. ASSOCIATED(scf_env%outer_scf%inv_jacobian)) &
ALLOCATE (scf_env%outer_scf%inv_jacobian(nvar, nvar))
inv_jacobian => scf_env%outer_scf%inv_jacobian
CALL invert_matrix(jacobian, inv_jacobian, inv_error)
DEALLOCATE (jacobian)
! Mark that an inverse Jacobian was just built and the next outer_loop_optimize should not perform
! a Broyden update of it
scf_control%outer_scf%broyden_update = .FALSE.
! Mark that the MD history has been used and should not be reused anymore on this energy evaluation
used_history = .TRUE.
END IF
END IF
IF (ihistory .GE. 2 .AND. .NOT. ASSOCIATED(scf_env%outer_scf%inv_jacobian)) THEN
! Next, try history from current SCF procedure
nvar = SIZE(scf_env%outer_scf%variables, 1)
IF (SIZE(scf_env%outer_scf%gradient, 2) .LT. 3) &
CALL cp_abort(__LOCATION__, &
"Keyword EXTRAPOLATION_ORDER in section OUTER_SCF must be greater than or equal "// &
"to 3 for optimizers that build the Jacobian from SCF history.")
ALLOCATE (jacobian(nvar, nvar))
DO i = 1, nvar
jacobian(i, i) = (scf_env%outer_scf%gradient(i, ihistory)-scf_env%outer_scf%gradient(i, ihistory-1))/ &
(scf_env%outer_scf%variables(i, ihistory)-scf_env%outer_scf%variables(i, ihistory-1))
END DO
IF (.NOT. ASSOCIATED(scf_env%outer_scf%inv_jacobian)) &
ALLOCATE (scf_env%outer_scf%inv_jacobian(nvar, nvar))
inv_jacobian => scf_env%outer_scf%inv_jacobian
CALL invert_matrix(jacobian, inv_jacobian, inv_error)
DEALLOCATE (jacobian)
scf_control%outer_scf%broyden_update = .FALSE.
ELSE
! No history => will fall back to SD optimizer in outer_loop_optimize
END IF
END IF
END IF
CALL timestop(handle)
END SUBROUTINE qs_calculate_inverse_jacobian
! **************************************************************************************************
!> \brief Perform backtracking line search to find the optimal step size for the CDFT constraint
!> optimizer. Assumes that the CDFT gradient function is a smooth function of the constraint
!> variables.
!> \param qs_env the qs_environment_type where to perform the line search
!> \par History
!> 02.2017 created [Nico Holmberg]
! **************************************************************************************************
SUBROUTINE qs_cdft_line_search(qs_env)
TYPE(qs_environment_type), POINTER :: qs_env
CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_cdft_line_search', &
routineP = moduleN//':'//routineN
INTEGER, PARAMETER :: max_linesearch = 5
CHARACTER(len=default_path_length) :: project_name
INTEGER :: handle, i, ispin, iter_count, lp, &
max_scf, nspins, nvar, output_unit
LOGICAL :: converged, do_linesearch, should_exit, &
should_stop
REAL(KIND=dp) :: alpha
REAL(KIND=dp), DIMENSION(:), POINTER :: energy
REAL(KIND=dp), DIMENSION(:, :), POINTER :: gradient, inv_jacobian
REAL(KIND=dp), EXTERNAL :: dnrm2
TYPE(cdft_control_type), POINTER :: cdft_control
TYPE(cp_logger_type), POINTER :: logger, tmp_logger
TYPE(cp_para_env_type), POINTER :: para_env
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: p_rmpv
TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao_kp
TYPE(dft_control_type), POINTER :: dft_control
TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mos, mos_stashed
TYPE(qs_energy_type), POINTER :: energy_qs
TYPE(qs_ks_env_type), POINTER :: ks_env
TYPE(qs_rho_type), POINTER :: rho
TYPE(qs_scf_env_type), POINTER :: scf_env
TYPE(scf_control_type), POINTER :: scf_control
CALL timeset(routineN, handle)
NULLIFY (energy, gradient, p_rmpv, rho_ao_kp, mos, rho, &
mos_stashed, ks_env, scf_env, scf_control, dft_control, &
cdft_control, inv_jacobian, para_env, &
tmp_logger, energy_qs)
logger => cp_get_default_logger()
CPASSERT(ASSOCIATED(qs_env))
CALL get_qs_env(qs_env, scf_env=scf_env, ks_env=ks_env, &
scf_control=scf_control, mos=mos, rho=rho, &
dft_control=dft_control, &
para_env=para_env, energy=energy_qs)
do_linesearch = .FALSE.
SELECT CASE (scf_control%outer_scf%optimizer)
CASE DEFAULT
do_linesearch = .FALSE.
CASE (outer_scf_optimizer_newton_ls)
do_linesearch = .TRUE.
CASE (outer_scf_optimizer_broyden)
SELECT CASE (scf_control%outer_scf%broyden_type)
CASE (broyden_type_1, broyden_type_2, broyden_type_1_explicit, broyden_type_2_explicit)
do_linesearch = .FALSE.
CASE (broyden_type_1_ls, broyden_type_1_explicit_ls, broyden_type_2_ls, broyden_type_2_explicit_ls)
cdft_control => dft_control%qs_control%cdft_control
IF (.NOT. ASSOCIATED(cdft_control)) &
CALL cp_abort(__LOCATION__, &
"Optimizers that perform a line search can"// &
" only be used together with a valid CDFT constraint")
IF (ASSOCIATED(scf_env%outer_scf%inv_jacobian)) &
do_linesearch = .TRUE.
END SELECT
END SELECT
IF (do_linesearch) THEN
cdft_control => dft_control%qs_control%cdft_control
IF (.NOT. ASSOCIATED(cdft_control)) &
CALL cp_abort(__LOCATION__, &
"Optimizers that perform a line search can"// &
" only be used together with a valid CDFT constraint")
CPASSERT(ASSOCIATED(scf_env%outer_scf%inv_jacobian))
alpha = 1.0_dp
iter_count = scf_env%outer_scf%iter_count
! Redirect output from line search procedure to a new file by creating a temporary logger
project_name = logger%iter_info%project_name
IF (para_env%mepos == para_env%source) THEN
lp = LEN_TRIM(project_name)
project_name(lp+1:LEN(project_name)) = "-LineSearch.out"
CALL open_file(file_name=project_name, file_status="UNKNOWN", &
file_action="WRITE", file_position="APPEND", &
unit_number=output_unit)
ELSE
output_unit = -1
END IF
CALL cp_logger_create(tmp_logger, &
para_env=para_env, &
default_global_unit_nr=output_unit, &
close_global_unit_on_dealloc=.FALSE.)
! Save last converged state so we can roll back to it (mo_coeff and some outer_loop variables)
nspins = dft_control%nspins
ALLOCATE (mos_stashed(nspins))
@ -1183,152 +1597,105 @@ CONTAINS
END DO
CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
p_rmpv => rho_ao_kp(:, 1)
cdft_control => dft_control%qs_control%cdft_control
IF (.NOT. ASSOCIATED(cdft_control)) &
CALL cp_abort(__LOCATION__, "Optimizers that need the explicit Jacobian can"// &
" only be used together with a valid CDFT constraint")
! Allocate work
nvar = SIZE(scf_env%outer_scf%variables, 1)
max_scf = scf_control%outer_scf%max_scf+1
iter_count = scf_env%outer_scf%iter_count
ALLOCATE (gradient(nvar, max_scf))
gradient = scf_env%outer_scf%gradient
ALLOCATE (energy(max_scf))
energy = scf_env%outer_scf%energy
ALLOCATE (jacobian(nvar, nvar))
jacobian = 0.0_dp
! Setup finite difference scheme
SELECT CASE (scf_control%outer_scf%jacobian_type)
CASE DEFAULT
CALL cp_abort(__LOCATION__, "Unknown Jacobian type: "// &
cp_to_string(scf_control%outer_scf%jacobian_type))
CASE (jacobian_fd1)
nwork = 0
pwork = 1
ALLOCATE (coeff(nwork:pwork), step_multiplier(nwork:pwork))
coeff(nwork) = -1.0_dp
coeff(pwork) = 1.0_dp
step_multiplier = 1.0_dp
dh = scf_control%outer_scf%jacobian_step
CASE (jacobian_fd1_backward)
nwork = -1
pwork = 0
ALLOCATE (coeff(nwork:pwork), step_multiplier(nwork:pwork))
coeff(nwork) = -1.0_dp
coeff(pwork) = 1.0_dp
step_multiplier = -1.0_dp
dh = scf_control%outer_scf%jacobian_step
CASE (jacobian_fd2)
nwork = 0
pwork = 2
ALLOCATE (coeff(nwork:pwork), step_multiplier(nwork:pwork))
coeff(0) = -3.0_dp/2.0_dp
coeff(1) = 2.0_dp
coeff(2) = -1.0_dp/2.0_dp
step_multiplier = 1.0_dp
step_multiplier(2) = 2.0_dp
dh = 2.0_dp*scf_control%outer_scf%jacobian_step
CASE (jacobian_fd2_backward)
nwork = -2
pwork = 0
ALLOCATE (coeff(nwork:pwork), step_multiplier(nwork:pwork))
coeff(0) = 3.0_dp/2.0_dp
coeff(-1) = -2.0_dp
coeff(-2) = 1.0_dp/2.0_dp
step_multiplier = -1.0_dp
step_multiplier(-2) = -2.0_dp
dh = 2.0_dp*scf_control%outer_scf%jacobian_step
CASE (jacobian_fd1_central)
nwork = -1
pwork = 1
ALLOCATE (coeff(nwork:pwork), step_multiplier(nwork:pwork))
coeff(nwork) = -1.0_dp/2.0_dp
coeff(pwork) = 1.0_dp/2.0_dp
step_multiplier = 1.0_dp
step_multiplier(nwork) = -1.0_dp
dh = 2.0_dp*scf_control%outer_scf%jacobian_step
END SELECT
DO i = 1, nvar
jacobian(i, :) = coeff(0)*scf_env%outer_scf%gradient(i, iter_count)
END DO
! Broyden optimizers: perform update of inv_jacobian if necessary
IF (scf_control%outer_scf%broyden_update) THEN
CALL outer_loop_optimize(scf_env, scf_control)
! Reset the variables and prevent a reupdate of inv_jacobian
scf_env%outer_scf%variables(:, iter_count+1) = 0
scf_control%outer_scf%broyden_update = .FALSE.
END IF
! Print some info
IF (output_unit > 0) THEN
WRITE (output_unit, FMT="(/,A)") &
" ============================== CDFT SCF LOOP =================================="
" ================================== LINE SEARCH STARTED ================================="
WRITE (output_unit, FMT="(A)") &
" Evaluating inverse Jacobian using finite differences"
SELECT CASE (scf_control%outer_scf%jacobian_type)
CASE (jacobian_fd1)
WRITE (output_unit, '(A)') " Type : First order forward difference"
CASE (jacobian_fd1_backward)
WRITE (output_unit, '(A)') " Type : First order backward difference"
CASE (jacobian_fd2)
WRITE (output_unit, '(A)') " Type : Second order forward difference"
CASE (jacobian_fd2_backward)
WRITE (output_unit, '(A)') " Type : Second order backward difference"
CASE (jacobian_fd1_central)
WRITE (output_unit, '(A)') " Type : First order central difference"
CASE DEFAULT
CALL cp_abort(__LOCATION__, "Unknown Jacobian type: "// &
cp_to_string(scf_control%outer_scf%jacobian_type))
END SELECT
WRITE (output_unit, '(A,ES12.4)') " Step size : ", scf_control%outer_scf%jacobian_step
" Evaluating optimal step size for optimizer"
WRITE (output_unit, '(A,I10,A,I10)') &
" Energy evaluation: ", cdft_control%ienergy, ", CDFT SCF iteration: ", iter_count
END IF
! Calculate the Jacobian by perturbing each Lagrangian and recalculating the energy self-consistently
DO i = 1, nvar
DO iwork = nwork, pwork
IF (iwork == 0) CYCLE
scf_env%outer_scf%variables(i, iter_count+1) = scf_env%outer_scf%variables(i, iter_count)+ &
step_multiplier(iwork)* &
scf_control%outer_scf%jacobian_step
CALL outer_loop_update_qs_env(qs_env, scf_env)
CALL qs_ks_did_change(ks_env, potential_changed=.TRUE.)
CALL outer_loop_switch(scf_env, scf_control, cdft_control, cdft2ot)
CALL scf_env_do_scf(scf_env=scf_env, scf_control=scf_control, qs_env=qs_env, &
converged=converged, should_stop=should_stop)
CALL outer_loop_switch(scf_env, scf_control, cdft_control, ot2cdft)
! Update (iter_count + 1) element of gradient
scf_env%outer_scf%iter_count = scf_env%outer_scf%iter_count+1
CALL outer_loop_gradient(qs_env, scf_env)
scf_env%outer_scf%iter_count = scf_env%outer_scf%iter_count-1
! Update Jacobian
DO j = 1, nvar
jacobian(j, i) = jacobian(j, i)+coeff(iwork)*scf_env%outer_scf%gradient(j, iter_count+1)
END DO
! Reset everything to last converged state
scf_env%outer_scf%variables(i, iter_count+1) = 0.0_dp
scf_env%outer_scf%gradient = gradient
scf_env%outer_scf%energy = energy
DO ispin = 1, nspins
CALL duplicate_mo_set(mos(ispin)%mo_set, mos_stashed(ispin)%mo_set)
IF (mos(ispin)%mo_set%use_mo_coeff_b) THEN
CALL copy_fm_to_dbcsr(mos(ispin)%mo_set%mo_coeff, &
mos(ispin)%mo_set%mo_coeff_b)
ENDIF
CALL calculate_density_matrix(mos(ispin)%mo_set, &
p_rmpv(ispin)%matrix)
END DO
CALL qs_rho_update_rho(rho, qs_env=qs_env)
CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.)
! Perform backtracking line search
CALL cp_add_default_logger(tmp_logger)
DO i = 1, max_linesearch
IF (output_unit > 0) THEN
WRITE (output_unit, FMT="(A)") " "
WRITE (output_unit, FMT="(A)") " #####################################"
WRITE (output_unit, '(A,I10,A)') &
" ### Line search step: ", i, " ###"
WRITE (output_unit, FMT="(A)") " #####################################"
END IF
inv_jacobian => scf_env%outer_scf%inv_jacobian
! Newton update of CDFT variables with a step size of alpha
scf_env%outer_scf%variables(:, iter_count+1) = scf_env%outer_scf%variables(:, iter_count)-alpha* &
MATMUL(inv_jacobian, scf_env%outer_scf%gradient(:, iter_count))
! With updated CDFT variables, perform SCF
CALL outer_loop_update_qs_env(qs_env, scf_env)
CALL qs_ks_did_change(ks_env, potential_changed=.TRUE.)
CALL outer_loop_switch(scf_env, scf_control, cdft_control, cdft2ot)
CALL scf_env_do_scf(scf_env=scf_env, scf_control=scf_control, qs_env=qs_env, &
converged=converged, should_stop=should_stop)
CALL outer_loop_switch(scf_env, scf_control, cdft_control, ot2cdft)
! Update (iter_count + 1) element of gradient and print constraint info
scf_env%outer_scf%iter_count = scf_env%outer_scf%iter_count+1
CALL outer_loop_gradient(qs_env, scf_env)
CALL qs_scf_cdft_info(output_unit, scf_control, scf_env, cdft_control, &
energy_qs, cdft_control%total_steps, &
should_stop=.FALSE., outer_loop_converged=.FALSE., cdft_loop=.FALSE.)
scf_env%outer_scf%iter_count = scf_env%outer_scf%iter_count-1
! Check if the L2 norm of the gradient decreased
inv_jacobian => scf_env%outer_scf%inv_jacobian
IF (dnrm2(nvar, scf_env%outer_scf%gradient(:, iter_count+1), 1) .LT. &
dnrm2(nvar, scf_env%outer_scf%gradient(:, iter_count), 1)) THEN
! Optimal step size found
should_exit = .TRUE.
ELSE
! Gradient increased => alpha is too large (if the gradient function is smooth)
should_exit = .FALSE.
! Update alpha using Armijo's scheme
alpha = alpha*0.5_dp
END IF
! Reset everything to last converged state
scf_env%outer_scf%variables(:, iter_count+1) = 0.0_dp
scf_env%outer_scf%gradient = gradient
scf_env%outer_scf%energy = energy
! If a suitable step size was found, do not reset the electronic structure so we can continue from it
IF (should_exit) EXIT
DO ispin = 1, nspins
CALL deallocate_mo_set(mos(ispin)%mo_set)
CALL duplicate_mo_set(mos(ispin)%mo_set, mos_stashed(ispin)%mo_set)
CALL calculate_density_matrix(mos(ispin)%mo_set, &
p_rmpv(ispin)%matrix)
END DO
CALL qs_rho_update_rho(rho, qs_env=qs_env)
CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.)
END DO
! Finalize and invert Jacobian
jacobian = jacobian/dh
IF (.NOT. ASSOCIATED(scf_env%outer_scf%inv_jacobian)) &
ALLOCATE (scf_env%outer_scf%inv_jacobian(nvar, nvar))
inv_jacobian => scf_env%outer_scf%inv_jacobian
CALL invert_matrix(jacobian, inv_jacobian, inv_error)
scf_control%outer_scf%newton_step = alpha
IF (.NOT. should_exit) THEN
CALL cp_warn(__LOCATION__, &
"Line search did not converge.")
scf_control%outer_scf%newton_step = 1.0_dp
END IF
CALL cp_rm_default_logger()
CALL cp_logger_release(tmp_logger)
! Release temporary storage
DO ispin = 1, nspins
CALL deallocate_mo_set(mos_stashed(ispin)%mo_set)
END DO
DEALLOCATE (mos_stashed, jacobian, gradient, energy, coeff, step_multiplier)
IF (output_unit > 0) &
DEALLOCATE (mos_stashed, gradient, energy)
IF (output_unit > 0) THEN
WRITE (output_unit, FMT="(/,A)") &
" ============================= JACOBIAN CALCULATED ============================="
" ================================== LINE SEARCH COMPLETE =================================="
END IF
END IF
CALL timestop(handle)
END SUBROUTINE outer_loop_calculate_inverse_jacobian
END SUBROUTINE qs_cdft_line_search
END MODULE qs_scf

22
src/qs_scf_initialization.F
View file

@ -308,8 +308,10 @@ CONTAINS
INTEGER :: nhistory, nvariables
LOGICAL :: do_kpoints
REAL(KIND=dp), DIMENSION(:, :), POINTER :: outer_scf_history
REAL(KIND=dp), DIMENSION(:, :), POINTER :: gradient_history, outer_scf_history, &
variable_history
NULLIFY (outer_scf_history, gradient_history, variable_history)
CALL get_qs_env(qs_env=qs_env, do_kpoints=do_kpoints)
! Check only one constraint section is active
IF (dft_control%qs_control%becke_restraint .AND. &
@ -343,7 +345,6 @@ CONTAINS
! Executed only on first call (OT settings active in scf_control)
! Save OT settings and constraint initial values in CDFT control
! Then switch to constraint outer_scf settings for proper initialization of history
NULLIFY (outer_scf_history)
IF (scf_control%outer_scf%have_scf) THEN
IF (scf_control%outer_scf%type == outer_scf_none) THEN
dft_control%qs_control%cdft_control%ot_control%have_scf = .TRUE.
@ -355,6 +356,12 @@ CONTAINS
dft_control%qs_control%cdft_control%ot_control%diis_buffer_length = scf_control%outer_scf%diis_buffer_length
dft_control%qs_control%cdft_control%ot_control%jacobian_type = scf_control%outer_scf%jacobian_type
dft_control%qs_control%cdft_control%ot_control%bisect_trust_count = scf_control%outer_scf%bisect_trust_count
dft_control%qs_control%cdft_control%ot_control%jacobian_step = scf_control%outer_scf%jacobian_step
dft_control%qs_control%cdft_control%ot_control%jacobian_freq = scf_control%outer_scf%jacobian_freq
dft_control%qs_control%cdft_control%ot_control%ijacobian = scf_control%outer_scf%ijacobian
dft_control%qs_control%cdft_control%ot_control%build_jacobian = scf_control%outer_scf%build_jacobian
dft_control%qs_control%cdft_control%ot_control%broyden_type = scf_control%outer_scf%broyden_type
dft_control%qs_control%cdft_control%ot_control%broyden_update = scf_control%outer_scf%broyden_update
nvariables = 1
! Constraint specific initializations
IF (dft_control%qs_control%becke_restraint) THEN
@ -369,14 +376,21 @@ CONTAINS
dft_control%qs_control%cdft_control%constraint_type = dft_control%qs_control%becke_control%constraint_type
dft_control%qs_control%cdft_control%combined_type = dft_control%qs_control%becke_control%combined_type
END IF
! In case constraint and ot extrapolation orders are different, make sure to use former
! In case constraint and OT extrapolation orders are different, make sure to use former
IF (scf_control%outer_scf%extrapolation_order /= &
dft_control%qs_control%cdft_control%constraint_control%extrapolation_order &
.OR. nvariables == 2) THEN
DEALLOCATE (qs_env%outer_scf_history)
DEALLOCATE (qs_env%gradient_history)
DEALLOCATE (qs_env%variable_history)
nhistory = dft_control%qs_control%cdft_control%constraint_control%extrapolation_order
ALLOCATE (outer_scf_history(nvariables, nhistory))
CALL set_qs_env(qs_env, outer_scf_history=outer_scf_history)
ALLOCATE (gradient_history(nvariables, 2))
gradient_history = 0.0_dp
ALLOCATE (variable_history(nvariables, 2))
variable_history = 0.0_dp
CALL set_qs_env(qs_env, outer_scf_history=outer_scf_history, &
gradient_history=gradient_history, variable_history=variable_history)
END IF
CALL outer_loop_switch(scf_env, scf_control, dft_control%qs_control%cdft_control, ot2cdft)
END IF

83
src/qs_scf_output.F
View file

@ -19,11 +19,15 @@ MODULE qs_scf_output
USE dbcsr_api, ONLY: dbcsr_p_type
USE input_constants, ONLY: &
becke_cavity_conf, becke_cutoff_element, becke_cutoff_global, becke_dynamic_conf, &
becke_mixed_conf, becke_none_conf, becke_static_conf, cdft_combined_acceptor, &
cdft_combined_all, cdft_combined_constraint, cdft_combined_donor, cdft_density_constraint, &
cdft_magnetization_constraint, outer_scf_becke_constraint, outer_scf_hirshfeld_constraint, &
outer_scf_optimizer_bisect, outer_scf_optimizer_broyden, outer_scf_optimizer_diis, &
outer_scf_optimizer_sd, radius_covalent, radius_single, radius_user, radius_vdw
becke_mixed_conf, becke_none_conf, becke_static_conf, broyden_type_1, &
broyden_type_1_explicit, broyden_type_1_explicit_ls, broyden_type_1_ls, broyden_type_2, &
broyden_type_2_explicit, broyden_type_2_explicit_ls, broyden_type_2_ls, &
cdft_combined_acceptor, cdft_combined_all, cdft_combined_constraint, cdft_combined_donor, &
cdft_density_constraint, cdft_magnetization_constraint, outer_scf_becke_constraint, &
outer_scf_hirshfeld_constraint, outer_scf_optimizer_bisect, outer_scf_optimizer_broyden, &
outer_scf_optimizer_diis, outer_scf_optimizer_newton, outer_scf_optimizer_newton_ls, &
outer_scf_optimizer_sd, outer_scf_optimizer_secant, radius_covalent, radius_single, &
radius_user, radius_vdw
USE input_section_types, ONLY: section_vals_get_subs_vals,&
section_vals_type,&
section_vals_val_get
@ -726,32 +730,30 @@ CONTAINS
SELECT CASE (cdft_control%constraint_type)
CASE (cdft_density_constraint)
WRITE (output_unit, '(A)') &
" Optimizing charge density constraint in an external SCF loop "
" Optimizing a charge density constraint in an external SCF loop "
CASE (cdft_magnetization_constraint)
WRITE (output_unit, '(A)') &
" Optimizing magnetization density constraint in an external SCF loop "
" Optimizing a magnetization density constraint in an external SCF loop "
CASE (cdft_combined_constraint)
WRITE (output_unit, '(A)') &
" Optimizing combined charge-spin constraint in an external SCF loop "
" Optimizing a combined charge-spin constraint in an external SCF loop "
SELECT CASE (cdft_control%combined_type)
CASE (cdft_combined_all)
WRITE (output_unit, '(A)') &
" spin constraint defined on all constraint atoms "
" the spin constraint is defined on all constraint atoms "
CASE (cdft_combined_acceptor)
WRITE (output_unit, '(A)') &
" spin constraint defined on acceptor constraint atoms "
" the spin constraint is defined on acceptor constraint atoms "
CASE (cdft_combined_donor)
WRITE (output_unit, '(A)') &
" spin constraint defined on donor constraint atoms "
" the spin constraint is defined on donor constraint atoms "
END SELECT
END SELECT
WRITE (output_unit, '(A)') " "
SELECT CASE (cdft_control%constraint_control%optimizer)
CASE (outer_scf_optimizer_sd)
WRITE (output_unit, '(A)') &
" Minimizer : SD : steepest descent"
CASE (outer_scf_optimizer_broyden)
WRITE (output_unit, '(A)') &
" Minimizer : Broyden : Broyden's method"
CASE (outer_scf_optimizer_diis)
WRITE (output_unit, '(A)') &
" Minimizer : DIIS : direct inversion"
@ -766,7 +768,60 @@ CONTAINS
WRITE (output_unit, '(A,I3)') &
" using a trust count of", &
cdft_control%constraint_control%bisect_trust_count
CASE (outer_scf_optimizer_broyden)
WRITE (output_unit, '(A)') " Minimizer : Broyden"
SELECT CASE (cdft_control%constraint_control%broyden_type)
CASE (broyden_type_1)
WRITE (output_unit, '(A)') " variant : 1st method"
CASE (broyden_type_1_explicit)
WRITE (output_unit, '(A)') " variant : 1st method with explicit initial Jacobian"
CASE (broyden_type_1_ls)
WRITE (output_unit, '(A)') " variant : 1st method with backtracking line search"
CASE (broyden_type_1_explicit_ls)
WRITE (output_unit, '(A)') &
" variant : 1st method with explicit initial Jacobian"
WRITE (output_unit, '(A)') &
" and backtracking line search"
CASE (broyden_type_2)
WRITE (output_unit, '(A)') " variant : 2nd method"
CASE (broyden_type_2_explicit)
WRITE (output_unit, '(A)') " variant : 2nd method with explicit initial Jacobian"
CASE (broyden_type_2_ls)
WRITE (output_unit, '(A)') " variant : 2nd method with backtracking line search"
CASE (broyden_type_2_explicit_ls)
WRITE (output_unit, '(A)') &
" variant : 2nd method with explicit initial Jacobian"
WRITE (output_unit, '(A)') &
" and backtracking line search"
END SELECT
WRITE (output_unit, '(A,I4,A)') &
" The Jacobian is restarted every ", &
cdft_control%constraint_control%jacobian_freq(2), " energy evaluation"
WRITE (output_unit, '(A,I4,A)') &
" or if the number of SCF iterations exceeds ", &
cdft_control%constraint_control%jacobian_freq(1)
CASE (outer_scf_optimizer_newton)
WRITE (output_unit, '(A)') " Minimizer : Newton"
WRITE (output_unit, '(A,I4,A)') &
" The Jacobian is restarted every ", &
cdft_control%constraint_control%jacobian_freq(2), " energy evaluation"
WRITE (output_unit, '(A,I4,A)') &
" or if the number of SCF iterations exceeds ", &
cdft_control%constraint_control%jacobian_freq(1)
CASE (outer_scf_optimizer_newton_ls)
WRITE (output_unit, '(A)') " Minimizer : Newton with backtracking line search"
WRITE (output_unit, '(A,I4,A)') &
" The Jacobian is restarted every ", &
cdft_control%constraint_control%jacobian_freq(2), " energy evaluation"
WRITE (output_unit, '(A,I4,A)') &
" or if the number of SCF iterations exceeds ", &
cdft_control%constraint_control%jacobian_freq(1)
CASE (outer_scf_optimizer_secant)
WRITE (output_unit, '(A)') " Minimizer : Secant"
CASE DEFAULT
CPABORT("")
END SELECT
WRITE (output_unit, '(A)') " "
SELECT CASE (cdft_control%type)
CASE (outer_scf_hirshfeld_constraint)
WRITE (output_unit, '(A)') " Type of constraint : Hirshfeld"

2
src/qs_scf_types.F
View file

@ -70,6 +70,7 @@ MODULE qs_scf_types
! **************************************************************************************************
TYPE qs_outer_scf_type
INTEGER :: iter_count
LOGICAL :: deallocate_jacobian
! these are the variable of outer loop.
! right now, we assume that they can be easily written as
! small arrays, but we might want to go the cp_fm_types
@ -173,6 +174,7 @@ CONTAINS
NULLIFY (scf_env%outer_scf%energy)
NULLIFY (scf_env%outer_scf%count)
NULLIFY (scf_env%outer_scf%inv_jacobian)
scf_env%outer_scf%deallocate_jacobian = .TRUE.
NULLIFY (scf_env%scf_work1)
NULLIFY (scf_env%scf_work2)
NULLIFY (scf_env%ortho)

88
src/scf_control_types.F
View file

@ -20,12 +20,15 @@ MODULE scf_control_types
cp_print_key_unit_nr
USE cp_units, ONLY: cp_unit_from_cp2k
USE input_constants, ONLY: &
atomic_guess, core_guess, diag_ot, direct_p_mix, general_roks, high_spin_roks, &
atomic_guess, broyden_type_1, broyden_type_1_explicit, broyden_type_1_explicit_ls, &
broyden_type_1_ls, broyden_type_2, broyden_type_2_explicit, broyden_type_2_explicit_ls, &
broyden_type_2_ls, core_guess, diag_ot, direct_p_mix, general_roks, high_spin_roks, &
ot_algo_taylor_or_diag, outer_scf_basis_center_opt, outer_scf_becke_constraint, &
outer_scf_ddapc_constraint, outer_scf_none, outer_scf_optimizer_bisect, &
outer_scf_optimizer_broyden, outer_scf_optimizer_diis, outer_scf_optimizer_newton, &
outer_scf_optimizer_none, outer_scf_optimizer_sd, outer_scf_s2_constraint, &
smear_energy_window, smear_fermi_dirac, smear_list
outer_scf_optimizer_newton_ls, outer_scf_optimizer_none, outer_scf_optimizer_sd, &
outer_scf_optimizer_secant, outer_scf_s2_constraint, smear_energy_window, &
smear_fermi_dirac, smear_list
USE input_cp2k_dft, ONLY: create_scf_section
USE input_enumeration_types, ONLY: enum_i2c,&
enumeration_type
@ -84,7 +87,13 @@ MODULE scf_control_types
!> \param build_jacobian logical which determines if the inverse Jacobian should be computed
!> \param step_size the optimizer step size
!> \param jacobian_step the step size for calculating the finite difference Jacobian
!> \param newton_step the step size used by Newton optimizer the with values between 0 and 1
!> \param jacobian_type the finite difference scheme to compute the Jacobian
!> \param broyden_type the variant of Broyden's method to use
!> \param jacobian_freq control parameters defining how often the Jacobian is built
!> \param ijacobian counter to track how many outer_scf iterations have passed since the last
!> Jacobian rebuild
!> \param broyden_update logical which determines if a Broyden update is needed
!> \par History
!> 09.2002 created [fawzi]
!> \author Fawzi Mohamed
@ -93,15 +102,20 @@ MODULE scf_control_types
TYPE outer_scf_control_type
LOGICAL :: have_scf
LOGICAL :: build_jacobian
LOGICAL :: broyden_update
INTEGER :: max_scf
REAL(KIND=dp) :: eps_scf, step_size, &
jacobian_step
jacobian_step, &
newton_step
INTEGER :: TYPE
INTEGER :: optimizer
INTEGER :: diis_buffer_length
INTEGER :: extrapolation_order
INTEGER :: bisect_trust_count
INTEGER :: jacobian_type
INTEGER :: broyden_type
INTEGER :: jacobian_freq(2)
INTEGER :: ijacobian
END TYPE outer_scf_control_type
TYPE smear_type
@ -268,8 +282,13 @@ CONTAINS
scf_control%outer_scf%optimizer = -1
scf_control%outer_scf%diis_buffer_length = -1
scf_control%outer_scf%jacobian_type = -1
scf_control%outer_scf%broyden_type = -1
scf_control%outer_scf%jacobian_freq = 1
scf_control%outer_scf%jacobian_step = 0.0_dp
scf_control%outer_scf%newton_step = 1.0_dp
scf_control%outer_scf%ijacobian = 0
scf_control%outer_scf%build_jacobian = .FALSE.
scf_control%outer_scf%broyden_update = .FALSE.
! Smearing of the MO occupations
@ -356,8 +375,8 @@ CONTAINS
routineP = moduleN//':'//routineN
INTEGER :: cholesky_flag, handle, ialgo
INTEGER, DIMENSION(:), POINTER :: added_mos
LOGICAL :: do_mixing
INTEGER, DIMENSION(:), POINTER :: added_mos, tmplist
LOGICAL :: do_mixing, exists
REAL(KIND=dp), DIMENSION(:), POINTER :: roks_parameter
TYPE(section_vals_type), POINTER :: mixing_section, outer_scf_section, &
scf_section, smear_section
@ -509,6 +528,23 @@ CONTAINS
i_val=scf_control%outer_scf%jacobian_type)
CALL section_vals_val_get(outer_scf_section, "JACOBIAN_STEP", &
r_val=scf_control%outer_scf%jacobian_step)
CALL section_vals_val_get(outer_scf_section, "BROYDEN_TYPE", &
i_val=scf_control%outer_scf%broyden_type)
CALL section_vals_val_get(outer_scf_section, "JACOBIAN_FREQ", explicit=exists)
IF (exists) THEN
CALL section_vals_val_get(outer_scf_section, "JACOBIAN_FREQ", &
i_vals=tmplist)
IF (SIZE(tmplist) /= 2) &
CALL cp_abort(__LOCATION__, &
"Keyword JACOBIAN_FREQ takes exactly two input values.")
IF (ANY(tmplist .LT. 0)) &
CALL cp_abort(__LOCATION__, &
"Keyword JACOBIAN_FREQ takes only positive values.")
IF (ALL(tmplist .EQ. 0)) &
CALL cp_abort(__LOCATION__, &
"Both values to keyword JACOBIAN_FREQ cannot be zero.")
scf_control%outer_scf%jacobian_freq(:) = tmplist(1:2)
END IF
END IF
smear_section => section_vals_get_subs_vals(scf_section, "SMEAR")
@ -755,8 +791,48 @@ CONTAINS
scf_control%outer_scf%diis_buffer_length
CASE (outer_scf_optimizer_broyden)
WRITE (output_unit, '(T25,A)') "Optimization with Broyden's method"
SELECT CASE (scf_control%outer_scf%broyden_type)
CASE (broyden_type_1)
WRITE (output_unit, '(A)') " variant : 1st method"
CASE (broyden_type_1_explicit)
WRITE (output_unit, '(A)') " variant : 1st method with explicit initial Jacobian"
CASE (broyden_type_1_ls)
WRITE (output_unit, '(A)') " variant : 1st method with backtracking line search"
CASE (broyden_type_1_explicit_ls)
WRITE (output_unit, '(A)') &
" variant : 1st method with explicit initial Jacobian"
WRITE (output_unit, '(A)') &
" and backtracking line search"
CASE (broyden_type_2)
WRITE (output_unit, '(A)') " variant : 2nd method"
CASE (broyden_type_2_explicit)
WRITE (output_unit, '(A)') " variant : 2nd method with explicit initial Jacobian"
CASE (broyden_type_2_ls)
WRITE (output_unit, '(A)') " variant : 2nd method with backtracking line search"
CASE (broyden_type_2_explicit_ls)
WRITE (output_unit, '(A)') &
" variant : 2nd method with explicit initial Jacobian"
WRITE (output_unit, '(A)') &
" and backtracking line search"
END SELECT
WRITE (output_unit, '(T25,A,I4,A)') &
"The Jacobian is restarted every ", scf_control%outer_scf%jacobian_freq(2), " energy evaluation"
WRITE (output_unit, '(T25,A,I4,A)') &
" or if the number of SCF iterations exceeds ", scf_control%outer_scf%jacobian_freq(1)
CASE (outer_scf_optimizer_newton)
WRITE (output_unit, '(T25,A)') "Optimization with Newton's method"
WRITE (output_unit, '(T25,A,I4,A)') &
"The Jacobian is restarted every ", scf_control%outer_scf%jacobian_freq(2), " energy evaluation"
WRITE (output_unit, '(T25,A,I4,A)') &
" or if the number of SCF iterations exceeds ", scf_control%outer_scf%jacobian_freq(1)
CASE (outer_scf_optimizer_newton_ls)
WRITE (output_unit, '(T25,A)') "Optimization with Newton's method using backtracking line search"
WRITE (output_unit, '(T25,A,I4,A)') &
"The Jacobian is restarted every ", scf_control%outer_scf%jacobian_freq(2), " energy evaluation"
WRITE (output_unit, '(T25,A,I4,A)') &
" or if the number of SCF iterations exceeds ", scf_control%outer_scf%jacobian_freq(1)
CASE (outer_scf_optimizer_secant)
WRITE (output_unit, '(T25,A)') "Optimization with the secant method"
CASE DEFAULT
CPABORT("")
END SELECT

48
tests/QMMM/QS/regtest-cdft/HeH-water-cdft-state-1.inp Normal file
View file

@ -0,0 +1,48 @@
@SET RESTART_WFN TRUE
@SET WFN_FILE HeH-water-noconstraint-1_0.wfn
@SET PROJECT_NAME HeH-water-cdft-state-1
@SET WRITE_WFN 1
@SET CHARGE 1
@SET WRITE_CUBE FALSE
@SET XYZFILE HeH-water.xyz
@SET BECKE_ACTIVE TRUE
@SET BECKE_FRAGMENT FALSE
@SET MAX_SCF 5
! He+ H
@SET BECKE_TARGET_1 0.0
@SET BECKE_STR_1 0.0
! He H+
@SET BECKE_TARGET_2 2.0
@SET BECKE_STR_2 0.0
@SET BECKE_GLOBAL_CUTOFF TRUE
@SET BECKE_CUTOFF_ELEMENT FALSE
@SET BECKE_ADJUST_SIZE FALSE
@SET BECKE_ATOMIC_CHARGES TRUE
@SET BECKE_CAVITY_CONFINE FALSE
@SET BECKE_CAVITY_SHAPE VDW
@SET BECKE_CAVITY_PRINT FALSE
@SET BECKE_SHOULD_SKIP FALSE
@SET BECKE_IN_MEMORY FALSE
&GLOBAL
PROJECT ${PROJECT_NAME}
RUN_TYPE ENERGY
PRINT_LEVEL MEDIUM
&END GLOBAL
&FORCE_EVAL
METHOD QMMM
@SET BECKE_STR ${BECKE_STR_1}
@SET BECKE_TARGET ${BECKE_TARGET_1}
@include subsys.inc
@include qmmm-common-params.inc
@include dft-common-params.inc
@include mm-common-params.inc
&END FORCE_EVAL

48
tests/QMMM/QS/regtest-cdft/HeH-water-cdft-state-2.inp Normal file
View file

@ -0,0 +1,48 @@
@SET RESTART_WFN TRUE
@SET WFN_FILE HeH-water-noconstraint-1_0.wfn
@SET PROJECT_NAME HeH-water-cdft-state-2
@SET WRITE_WFN 1
@SET CHARGE 1
@SET WRITE_CUBE FALSE
@SET XYZFILE HeH-water.xyz
@SET BECKE_ACTIVE TRUE
@SET BECKE_FRAGMENT FALSE
@SET MAX_SCF 5
! He+ H
@SET BECKE_TARGET_1 0.0
@SET BECKE_STR_1 0.0
! He H+
@SET BECKE_TARGET_2 2.0
@SET BECKE_STR_2 0.0
@SET BECKE_GLOBAL_CUTOFF TRUE
@SET BECKE_CUTOFF_ELEMENT FALSE
@SET BECKE_ADJUST_SIZE FALSE
@SET BECKE_ATOMIC_CHARGES TRUE
@SET BECKE_CAVITY_CONFINE FALSE
@SET BECKE_CAVITY_SHAPE VDW
@SET BECKE_CAVITY_PRINT FALSE
@SET BECKE_SHOULD_SKIP FALSE
@SET BECKE_IN_MEMORY FALSE
&GLOBAL
PROJECT ${PROJECT_NAME}
RUN_TYPE ENERGY
PRINT_LEVEL MEDIUM
&END GLOBAL
&FORCE_EVAL
METHOD QMMM
@SET BECKE_STR ${BECKE_STR_2}
@SET BECKE_TARGET ${BECKE_TARGET_2}
@include subsys.inc
@include qmmm-common-params.inc
@include dft-common-params.inc
@include mm-common-params.inc
&END FORCE_EVAL

97
tests/QMMM/QS/regtest-cdft/HeH-water-mixed-cdft-1.inp Normal file
View file

@ -0,0 +1,97 @@
@SET RESTART_WFN TRUE
@SET WFN_FILE_1 HeH-water-cdft-state-1-1_0.wfn
@SET WFN_FILE_2 HeH-water-cdft-state-2-1_0.wfn
@SET PROJECT_NAME HeH-water-mixed-cdft-1
@SET NAME ${PROJECT_NAME}
@SET WRITE_WFN 0
@SET CHARGE 1
@SET WRITE_CUBE FALSE
@SET XYZFILE HeH-water.xyz
@SET BECKE_ACTIVE TRUE
@SET BECKE_FRAGMENT FALSE
@SET MAX_SCF 5
! He+ H
@SET BECKE_TARGET_1 0.0
@SET BECKE_STR_1 0.575328339408
! He H+
@SET BECKE_TARGET_2 2.0
@SET BECKE_STR_2 -1.885819092012
@SET BECKE_GLOBAL_CUTOFF TRUE
@SET BECKE_CUTOFF_ELEMENT FALSE
@SET BECKE_ADJUST_SIZE FALSE
@SET BECKE_ATOMIC_CHARGES FALSE
@SET BECKE_CAVITY_CONFINE TRUE
@SET BECKE_CAVITY_SHAPE VDW
@SET BECKE_CAVITY_PRINT FALSE
@SET BECKE_SHOULD_SKIP FALSE
@SET BECKE_IN_MEMORY TRUE
@SET USE_DLB FALSE
@SET COMPUTE_METRIC FALSE
@SET WFN_OVERLAP_METHOD FALSE
@SET NGROUPS 2
&GLOBAL
PROJECT ${PROJECT_NAME}
RUN_TYPE ENERGY_FORCE
PRINT_LEVEL MEDIUM
&END GLOBAL
&MULTIPLE_FORCE_EVALS
FORCE_EVAL_ORDER 2 3
MULTIPLE_SUBSYS F
&END
&FORCE_EVAL
METHOD MIXED
&MIXED
MIXING_TYPE LINEAR_COMBINATION
NGROUPS ${NGROUPS}
&LINEAR
LAMBDA 1.0
MIXED_CDFT TRUE
MIXED_CDFT_COUPLING 1
@IF ( ${USE_DLB} == TRUE )
! Defaults to FALSE
MIXED_CDFT_DLB TRUE
LOAD_SCALE 2
MORE_WORK 0
VERY_OVERLOADED 15
@ENDIF
&END LINEAR
&PRINT
&PROGRAM_RUN_INFO
&END
&END PRINT
&END MIXED
@include subsys.inc
&END FORCE_EVAL
&FORCE_EVAL
METHOD QMMM
@SET BECKE_STR ${BECKE_STR_1}
@SET BECKE_TARGET ${BECKE_TARGET_1}
@SET PROJECT_NAME ${NAME}-state-1
@SET WFN_FILE ${WFN_FILE_1}
@include qmmm-common-params.inc
@include dft-common-params.inc
@include mm-common-params.inc
&END FORCE_EVAL
&FORCE_EVAL
METHOD QMMM
@SET BECKE_STR ${BECKE_STR_2}
@SET BECKE_TARGET ${BECKE_TARGET_2}
@SET PROJECT_NAME ${NAME}-state-2
@SET WFN_FILE ${WFN_FILE_2}
@include qmmm-common-params.inc
@include dft-common-params.inc
@include mm-common-params.inc
&END FORCE_EVAL

97
tests/QMMM/QS/regtest-cdft/HeH-water-mixed-cdft-2.inp Normal file
View file

@ -0,0 +1,97 @@
@SET RESTART_WFN TRUE
@SET WFN_FILE_1 HeH-water-cdft-state-1-1_0.wfn
@SET WFN_FILE_2 HeH-water-cdft-state-2-1_0.wfn
@SET PROJECT_NAME HeH-water-mixed-cdft-2
@SET NAME ${PROJECT_NAME}
@SET WRITE_WFN 0
@SET CHARGE 1
@SET WRITE_CUBE FALSE
@SET XYZFILE HeH-water.xyz
@SET BECKE_ACTIVE TRUE
@SET BECKE_FRAGMENT FALSE
@SET MAX_SCF 5
! He+ H
@SET BECKE_TARGET_1 0.0
@SET BECKE_STR_1 0.575328339408
! He H+
@SET BECKE_TARGET_2 2.0
@SET BECKE_STR_2 -1.885819092012
@SET BECKE_GLOBAL_CUTOFF TRUE
@SET BECKE_CUTOFF_ELEMENT FALSE
@SET BECKE_ADJUST_SIZE FALSE
@SET BECKE_ATOMIC_CHARGES FALSE
@SET BECKE_CAVITY_CONFINE TRUE
@SET BECKE_CAVITY_SHAPE VDW
@SET BECKE_CAVITY_PRINT FALSE
@SET BECKE_SHOULD_SKIP FALSE
@SET BECKE_IN_MEMORY TRUE
@SET USE_DLB TRUE
@SET COMPUTE_METRIC FALSE
@SET WFN_OVERLAP_METHOD FALSE
@SET NGROUPS 1
&GLOBAL
PROJECT ${PROJECT_NAME}
RUN_TYPE ENERGY_FORCE
PRINT_LEVEL MEDIUM
&END GLOBAL
&MULTIPLE_FORCE_EVALS
FORCE_EVAL_ORDER 2 3
MULTIPLE_SUBSYS F
&END
&FORCE_EVAL
METHOD MIXED
&MIXED
MIXING_TYPE LINEAR_COMBINATION
NGROUPS ${NGROUPS}
&LINEAR
LAMBDA 1.0
MIXED_CDFT TRUE
MIXED_CDFT_COUPLING 1
@IF ( ${USE_DLB} == TRUE )
! Defaults to FALSE
MIXED_CDFT_DLB TRUE
LOAD_SCALE 2
MORE_WORK 0
VERY_OVERLOADED 15
@ENDIF
&END LINEAR
&PRINT
&PROGRAM_RUN_INFO
&END
&END PRINT
&END MIXED
@include subsys.inc
&END FORCE_EVAL
&FORCE_EVAL
METHOD QMMM
@SET BECKE_STR ${BECKE_STR_1}
@SET BECKE_TARGET ${BECKE_TARGET_1}
@SET PROJECT_NAME ${NAME}-state-1
@SET WFN_FILE ${WFN_FILE_1}
@include qmmm-common-params.inc
@include dft-common-params.inc
@include mm-common-params.inc
&END FORCE_EVAL
&FORCE_EVAL
METHOD QMMM
@SET BECKE_STR ${BECKE_STR_2}
@SET BECKE_TARGET ${BECKE_TARGET_2}
@SET PROJECT_NAME ${NAME}-state-2
@SET WFN_FILE ${WFN_FILE_2}
@include qmmm-common-params.inc
@include dft-common-params.inc
@include mm-common-params.inc
&END FORCE_EVAL

109
tests/QMMM/QS/regtest-cdft/HeH-water-mixed-cdft-3.inp Normal file
View file

@ -0,0 +1,109 @@
@SET RESTART_WFN TRUE
@SET WFN_FILE_1 HeH-water-cdft-state-1-1_0.wfn
@SET WFN_FILE_2 HeH-water-cdft-state-2-1_0.wfn
@SET PROJECT_NAME HeH-water-mixed-cdft-3
@SET NAME ${PROJECT_NAME}
@SET WRITE_WFN 0
@SET CHARGE 1
@SET WRITE_CUBE FALSE
@SET XYZFILE HeH-water.xyz
@SET BECKE_ACTIVE TRUE
@SET BECKE_FRAGMENT FALSE
@SET MAX_SCF 5
! He+ H
@SET BECKE_TARGET_1 0.0
@SET BECKE_STR_1 0.575328339408
! He H+
@SET BECKE_TARGET_2 2.0
@SET BECKE_STR_2 -1.885819092012
@SET BECKE_GLOBAL_CUTOFF TRUE
@SET BECKE_CUTOFF_ELEMENT FALSE
@SET BECKE_ADJUST_SIZE FALSE
@SET BECKE_ATOMIC_CHARGES FALSE
@SET BECKE_CAVITY_CONFINE TRUE
@SET BECKE_CAVITY_SHAPE VDW
@SET BECKE_CAVITY_PRINT FALSE
@SET BECKE_SHOULD_SKIP FALSE
@SET BECKE_IN_MEMORY TRUE
@SET USE_DLB FALSE
@SET COMPUTE_METRIC FALSE
@SET WFN_OVERLAP_METHOD FALSE
@SET NGROUPS 2
&GLOBAL
PROJECT ${PROJECT_NAME}
RUN_TYPE MD
PRINT_LEVEL MEDIUM
&END GLOBAL
&MOTION
&MD
ENSEMBLE NVT
STEPS 2
TIMESTEP 0.5
TEMPERATURE 300
&THERMOSTAT
TYPE CSVR
&END THERMOSTAT
&END MD
&END MOTION
&MULTIPLE_FORCE_EVALS
FORCE_EVAL_ORDER 2 3
MULTIPLE_SUBSYS F
&END
&FORCE_EVAL
METHOD MIXED
&MIXED
MIXING_TYPE LINEAR_COMBINATION
NGROUPS ${NGROUPS}
&LINEAR
LAMBDA 1.0
MIXED_CDFT TRUE
MIXED_CDFT_COUPLING 1
@IF ( ${USE_DLB} == TRUE )
! Defaults to FALSE
MIXED_CDFT_DLB TRUE
LOAD_SCALE 2
MORE_WORK 0
VERY_OVERLOADED 15
@ENDIF
&END LINEAR
&PRINT
&PROGRAM_RUN_INFO
&END
&END PRINT
&END MIXED
@include subsys.inc
&END FORCE_EVAL
&FORCE_EVAL
METHOD QMMM
@SET BECKE_STR ${BECKE_STR_1}
@SET BECKE_TARGET ${BECKE_TARGET_1}
@SET PROJECT_NAME ${NAME}-state-1
@SET WFN_FILE ${WFN_FILE_1}
@include qmmm-common-params.inc
@include dft-common-params.inc
@include mm-common-params.inc
&END FORCE_EVAL
&FORCE_EVAL
METHOD QMMM
@SET BECKE_STR ${BECKE_STR_2}
@SET BECKE_TARGET ${BECKE_TARGET_2}
@SET PROJECT_NAME ${NAME}-state-2
@SET WFN_FILE ${WFN_FILE_2}
@include qmmm-common-params.inc
@include dft-common-params.inc
@include mm-common-params.inc
&END FORCE_EVAL

109
tests/QMMM/QS/regtest-cdft/HeH-water-mixed-cdft-4.inp Normal file
View file

@ -0,0 +1,109 @@
@SET RESTART_WFN TRUE
@SET WFN_FILE_1 HeH-water-cdft-state-1-1_0.wfn
@SET WFN_FILE_2 HeH-water-cdft-state-2-1_0.wfn
@SET PROJECT_NAME HeH-water-mixed-cdft-4
@SET NAME ${PROJECT_NAME}
@SET WRITE_WFN 0
@SET CHARGE 1
@SET WRITE_CUBE FALSE
@SET XYZFILE HeH-water.xyz
@SET BECKE_ACTIVE TRUE
@SET BECKE_FRAGMENT FALSE
@SET MAX_SCF 5
! He+ H
@SET BECKE_TARGET_1 0.0
@SET BECKE_STR_1 0.575328339408
! He H+
@SET BECKE_TARGET_2 2.0
@SET BECKE_STR_2 -1.885819092012
@SET BECKE_GLOBAL_CUTOFF TRUE
@SET BECKE_CUTOFF_ELEMENT FALSE
@SET BECKE_ADJUST_SIZE FALSE
@SET BECKE_ATOMIC_CHARGES FALSE
@SET BECKE_CAVITY_CONFINE TRUE
@SET BECKE_CAVITY_SHAPE VDW
@SET BECKE_CAVITY_PRINT FALSE
@SET BECKE_SHOULD_SKIP FALSE
@SET BECKE_IN_MEMORY TRUE
@SET USE_DLB TRUE
@SET COMPUTE_METRIC FALSE
@SET WFN_OVERLAP_METHOD FALSE
@SET NGROUPS 1
&GLOBAL
PROJECT ${PROJECT_NAME}
RUN_TYPE MD
PRINT_LEVEL MEDIUM
&END GLOBAL
&MOTION
&MD
ENSEMBLE NVT
STEPS 2
TIMESTEP 0.5
TEMPERATURE 300
&THERMOSTAT
TYPE CSVR
&END THERMOSTAT
&END MD
&END MOTION
&MULTIPLE_FORCE_EVALS
FORCE_EVAL_ORDER 2 3
MULTIPLE_SUBSYS F
&END
&FORCE_EVAL
METHOD MIXED
&MIXED
MIXING_TYPE LINEAR_COMBINATION
NGROUPS ${NGROUPS}
&LINEAR
LAMBDA 1.0
MIXED_CDFT TRUE
MIXED_CDFT_COUPLING 1
@IF ( ${USE_DLB} == TRUE )
! Defaults to FALSE
MIXED_CDFT_DLB TRUE
LOAD_SCALE 2
MORE_WORK 0
VERY_OVERLOADED 15
@ENDIF
&END LINEAR
&PRINT
&PROGRAM_RUN_INFO
&END
&END PRINT
&END MIXED
@include subsys.inc
&END FORCE_EVAL
&FORCE_EVAL
METHOD QMMM
@SET BECKE_STR ${BECKE_STR_1}
@SET BECKE_TARGET ${BECKE_TARGET_1}
@SET PROJECT_NAME ${NAME}-state-1
@SET WFN_FILE ${WFN_FILE_1}
@include qmmm-common-params.inc
@include dft-common-params.inc
@include mm-common-params.inc
&END FORCE_EVAL
&FORCE_EVAL
METHOD QMMM
@SET BECKE_STR ${BECKE_STR_2}
@SET BECKE_TARGET ${BECKE_TARGET_2}
@SET PROJECT_NAME ${NAME}-state-2
@SET WFN_FILE ${WFN_FILE_2}
@include qmmm-common-params.inc
@include dft-common-params.inc
@include mm-common-params.inc
&END FORCE_EVAL

36
tests/QMMM/QS/regtest-cdft/HeH-water-noconstraint-dump-1.psf Normal file
View file

@ -0,0 +1,36 @@
PSF EXT
1 !NTITLE
CP2K generated DUMP of connectivity
5 !NATOM
1 MOL1 1 R1 OW OW -0.834000 15.999 0
2 MOL1 1 R1 HW HW 0.417000 1.008 0
3 MOL1 1 R1 HW HW 0.417000 1.008 0
4 MOL2 1 R2 He He -99.000000 4.003 0
5 MOL3 1 R3 H H -99.000000 1.008 0
2 !NBOND
1 2 1 3
1 !NTHETA
2 1 3
0 !NPHI
0 !NIMPHI
0 !NDON
0 !NACC
0 !NNB

46
tests/QMMM/QS/regtest-cdft/HeH-water-noconstraint.inp Normal file
View file

@ -0,0 +1,46 @@
@SET RESTART_WFN FALSE
@SET WFN_FILE HeH-water-noconstraint-1_0.wfn
@SET PROJECT_NAME HeH-water-noconstraint
@SET WRITE_WFN 1
@SET CHARGE 1
@SET WRITE_CUBE FALSE
@SET XYZFILE HeH-water.xyz
@SET BECKE_ACTIVE FALSE
@SET BECKE_FRAGMENT FALSE
@SET MAX_SCF 0
! He+ H
@SET BECKE_TARGET_1 0.0
@SET BECKE_STR_1 0.0
! He H+
@SET BECKE_TARGET_2 2.0
@SET BECKE_STR_2 0.0
@SET BECKE_GLOBAL_CUTOFF TRUE
@SET BECKE_CUTOFF_ELEMENT FALSE
@SET BECKE_ADJUST_SIZE FALSE
@SET BECKE_ATOMIC_CHARGES FALSE
@SET BECKE_CAVITY_CONFINE FALSE
@SET BECKE_CAVITY_SHAPE VDW
@SET BECKE_CAVITY_PRINT FALSE
@SET BECKE_SHOULD_SKIP FALSE
@SET BECKE_IN_MEMORY FALSE
&GLOBAL
PROJECT ${PROJECT_NAME}
RUN_TYPE ENERGY
PRINT_LEVEL MEDIUM
&END GLOBAL
&FORCE_EVAL
METHOD QMMM
@include subsys.inc
@include qmmm-common-params.inc
@include dft-common-params.inc
@include mm-common-params.inc
&END FORCE_EVAL

7
tests/QMMM/QS/regtest-cdft/HeH-water.xyz Normal file
View file

@ -0,0 +1,7 @@
5
OW 0.000000 0.000000 0.000000
HW 0.000000 0.950000 0.000000
HW 0.950000 0.000000 0.000000
He 7.500000 7.500000 7.100000
H 7.500000 7.500000 7.900000

17
tests/QMMM/QS/regtest-cdft/TEST_FILES Normal file
View file

@ -0,0 +1,17 @@
# 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-water-noconstraint.inp 1 1e-11 -2.95332848940723
# The MM part seems to introduce noise that is difficult to account for
# So just running the test succesfully is sufficient
HeH-water-cdft-state-1.inp 0
HeH-water-cdft-state-2.inp 0
# These tests compute the electronic coupling in parallel/serial modes
# Tests the charge transfer energy which is not prone to noise
HeH-water-mixed-cdft-1.inp 75 5e-10 0.318940758387
HeH-water-mixed-cdft-2.inp 75 5e-10 0.318940758387
# MD tests in parallel/serial modes
HeH-water-mixed-cdft-3.inp 11 5e-10 -2.564762048904492
HeH-water-mixed-cdft-4.inp 11 5e-10 -2.564762048791228
#EOF

3
tests/QMMM/QS/regtest-cdft/TEST_FILES_RESET Normal file
View file

@ -0,0 +1,3 @@
#
# add files to be reset here
#

88
tests/QMMM/QS/regtest-cdft/becke_qs.inc Normal file
View file

@ -0,0 +1,88 @@
&QS
METHOD GPW
EPS_DEFAULT 1.0E-12
MAP_CONSISTENT
EXTRAPOLATION ASPC
EXTRAPOLATION_ORDER 3
&CDFT
TYPE_OF_CONSTRAINT BECKE
&OUTER_SCF ON
EPS_SCF 3.0e-1
TYPE BECKE_CONSTRAINT
OPTIMIZER NEWTON
BROYDEN_TYPE BT1
JACOBIAN_STEP 1.0E-4
JACOBIAN_TYPE FD1
JACOBIAN_FREQ 1 1
BISECT_TRUST_COUNT 8
EXTRAPOLATION_ORDER 2
MAX_SCF ${MAX_SCF}
STEP_SIZE -1.0
&END
&END CDFT
&BECKE_RESTRAINT
@IF ( ${BECKE_ADJUST_SIZE} == TRUE )
! Defaults to false
ADJUST_SIZE TRUE
ATOMIC_RADII 0.460 0.320
@ENDIF
@IF ( ${BECKE_ATOMIC_CHARGES} == TRUE )
! Defaults to false
ATOMIC_CHARGES TRUE
@ENDIF
STRENGTH ${BECKE_STR}
! The constraint target: sum_i coeff_i * N_i
! where N_i is the number of VALENCE electrons on i
TARGET ${BECKE_TARGET}
@IF ( ${BECKE_CUTOFF_ELEMENT} == TRUE )
CUTOFF_TYPE ELEMENT
ELEMENT_CUTOFF 2.0 2.0
@ENDIF
@IF ( ${BECKE_GLOBAL_CUTOFF} == TRUE )
CUTOFF_TYPE GLOBAL
GLOBAL_CUTOFF 2.0
@ENDIF
@IF ( ${BECKE_IN_MEMORY} == TRUE )
! Defaults to false
IN_MEMORY TRUE
@ENDIF
@IF ( ${BECKE_CAVITY_CONFINE} == TRUE )
! Defaults to NONE
CONFINE CAVITY
EPS_CAVITY 1.0E-6
CAVITY_SHAPE ${BECKE_CAVITY_SHAPE}
! For shape single
CAVITY_RADIUS 1.3
CAVITY_USE_BOHR FALSE
@ENDIF
@IF ( ${BECKE_SHOULD_SKIP} == TRUE )
! Defaults to false
SHOULD_SKIP TRUE
@ENDIF
@IF ( ${BECKE_CAVITY_PRINT} == TRUE )
! Defaults to false
CAVITY_PRINT TRUE
@ENDIF
@IF ( ${BECKE_FRAGMENT} == TRUE )
! Need absolute weight to use fragment densities
ATOMS 1
COEFF 1
FRAGMENT_DENSITIES
FRAGMENT_A_FILE_NAME He+-noconstraint-ELECTRON_DENSITY-1_0.cube
FRAGMENT_B_FILE_NAME H-noconstraint-ELECTRON_DENSITY-1_0.cube
@ENDIF
@IF ( ${BECKE_FRAGMENT} == FALSE )
ATOMS 1..2
COEFF 1 -1
@ENDIF
CONSTRAINT_TYPE TOTAL
&PROGRAM_RUN_INFO ON
&EACH
QS_SCF 1
&END EACH
COMMON_ITERATION_LEVELS 2
ADD_LAST NUMERIC
FILENAME ./${PROJECT_NAME}
&END PROGRAM_RUN_INFO
&END BECKE_RESTRAINT
&END QS

83
tests/QMMM/QS/regtest-cdft/dft-common-params.inc Normal file
View file

@ -0,0 +1,83 @@
&DFT
@IF ( ${BECKE_ACTIVE} == TRUE )
@include becke_qs.inc
@ENDIF
@IF ( ${BECKE_ACTIVE} == FALSE )
&QS
METHOD GPW
EPS_DEFAULT 1.0E-12
MAP_CONSISTENT
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
COMMENSURATE
&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
&END
&END

71
tests/QMMM/QS/regtest-cdft/mm-common-params.inc Normal file
View file

@ -0,0 +1,71 @@
&MM
&FORCEFIELD
IGNORE_MISSING_CRITICAL_PARAMS
&SPLINE
RCUT_NB 1.2500000000000000E+01
EMAX_SPLINE 1.0000000000000000E+02
&END SPLINE
&NONBONDED
&LENNARD-JONES
ATOMS OW OW
EPSILON [kjmol] 0.63639
SIGMA [nm] 0.315057422683
&END LENNARD-JONES
&LENNARD-JONES
ATOMS OW HW
EPSILON [kjmol] 0.3499745204439889
SIGMA [nm] 0.17752938756375
&END LENNARD-JONES
&LENNARD-JONES
ATOMS HW HW
EPSILON [kjmol] 0.192464
SIGMA [nm] 0.0400013524445
&END LENNARD-JONES
&END NONBONDED
&CHARGE
ATOM OW
CHARGE -0.834
&END CHARGE
&CHARGE
ATOM HW
CHARGE 0.417
&END CHARGE
&BOND
ATOMS OW HW
KIND HARMONIC
K [kjmol*nm^-2] 376560.0
R0 [nm] 0.09572
&END BOND
&BOND
! Dummy
ATOMS HW HW
KIND HARMONIC
K 0.0000000000000000E+00
R0 0.09572
&END BOND
&BEND
ATOMS HW OW HW
KIND HARMONIC
K [kjmol*rad^-2] 460.24
THETA0 [deg] 104.52
&END BEND
&BEND
! Dummy
ATOMS OW HW HW
KIND HARMONIC
K 0.0000000000000000E+00
THETA0 1.9111355309337910E+00
&END BEND
&END FORCEFIELD
&POISSON
&EWALD
EWALD_TYPE spme
ALPHA 0.44
GMAX 30
&END EWALD
&END POISSON
&PRINT
&FF_INFO ON
&END FF_INFO
&END PRINT
&END MM

50
tests/QMMM/QS/regtest-cdft/qmmm-common-params.inc Normal file
View file

@ -0,0 +1,50 @@
&QMMM
USE_GEEP_LIB 4
CENTER EVERY_STEP
CENTER_GRID TRUE
NOCOMPATIBILITY
&PRINT
&MM_POTENTIAL
&END MM_POTENTIAL
&END
&PERIODIC
&MULTIPOLE
RCUT [angstrom] 3.5
ANALYTICAL_GTERM T
&END MULTIPOLE
&POISSON
PERIODIC XYZ
&END POISSON
&END
&CELL
ABC 8. 8. 8.
PERIODIC XYZ
&END CELL
ECOUPL GAUSS
&MM_KIND HW
RADIUS 0.440
&END MM_KIND
&MM_KIND OW
RADIUS 1.20
&END MM_KIND
&QM_KIND He
MM_INDEX 4
&END QM_KIND
&QM_KIND H
MM_INDEX 5
&END QM_KIND
&FORCEFIELD
&NONBONDED
&LENNARD-JONES
ATOMS HW O
EPSILON [kjmol] 0.4112256910262295
SIGMA [nm] 0.16799667622224998
&END LENNARD-JONES
&LENNARD-JONES
ATOMS HW H
EPSILON [kjmol] 0.109904688889965
SIGMA [nm] 0.14998267622224998
&END LENNARD-JONES
&END NONBONDED
&END
&END QMMM

26
tests/QMMM/QS/regtest-cdft/subsys.inc Normal file
View file

@ -0,0 +1,26 @@
&SUBSYS
&CELL
ABC 15. 15. 15.
PERIODIC XYZ
&END CELL
&TOPOLOGY
COORD_FILE_NAME ${XYZFILE}
COORD_FILE_FORMAT XYZ
CONN_FILE_FORMAT PSF
CONN_FILE_NAME HeH-water-noconstraint-dump-1.psf
&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
&KIND OW
ELEMENT O
&END KIND
&KIND HW
ELEMENT H
&END KIND
&END SUBSYS

4
tests/QS/regtest-cdft-3/HeH-mixed-cdft-1.inp
View file

@ -37,6 +37,8 @@
@SET COMPUTE_METRIC FALSE
@SET WFN_OVERLAP_METHOD FALSE
@SET NGROUPS 2
&GLOBAL
PROJECT ${PROJECT_NAME}
RUN_TYPE ENERGY
@ -51,7 +53,7 @@
METHOD MIXED
&MIXED
MIXING_TYPE LINEAR_COMBINATION
NGROUPS 2
NGROUPS ${NGROUPS}
&LINEAR
LAMBDA 1.0
MIXED_CDFT TRUE

116
tests/QS/regtest-cdft-3/HeH-mixed-cdft-10.inp Normal file
View file

@ -0,0 +1,116 @@
@SET RESTART_WFN TRUE
@SET WFN_FILE_1 HeH-cdft-state-1-1_0.wfn
@SET WFN_FILE_2 HeH-cdft-state-2-1_0.wfn
@SET PROJECT_NAME HeH-mixed-cdft-10
@SET NAME ${PROJECT_NAME}
@SET WRITE_WFN 0
@SET CHARGE 1
@SET WRITE_CUBE FALSE
@SET XYZFILE HeH.xyz
@SET BECKE_ACTIVE TRUE
@SET BECKE_FRAGMENT FALSE
@SET MAX_SCF 5
! He+ H
@SET BECKE_TARGET_1 0.0
@SET BECKE_STR_1 0.707711420058
! He H+
@SET BECKE_TARGET_2 2.0
@SET BECKE_STR_2 -1.008285468010
@SET BECKE_GLOBAL_CUTOFF FALSE
@SET BECKE_CUTOFF_ELEMENT TRUE
@SET BECKE_ADJUST_SIZE TRUE
@SET BECKE_ATOMIC_CHARGES FALSE
@SET BECKE_CAVITY_CONFINE TRUE
@SET BECKE_CAVITY_SHAPE VDW
@SET BECKE_CAVITY_PRINT FALSE
@SET BECKE_SHOULD_SKIP TRUE
@SET BECKE_IN_MEMORY TRUE
@SET USE_DLB FALSE
@SET COMPUTE_METRIC FALSE
@SET WFN_OVERLAP_METHOD FALSE
@SET NGROUPS 1
&GLOBAL
PROJECT ${PROJECT_NAME}
RUN_TYPE MD
PRINT_LEVEL MEDIUM
&TIMINGS
THRESHOLD 0.001
&END
&END GLOBAL
&MOTION
&MD
ENSEMBLE NVT
STEPS 2
TIMESTEP 0.5
TEMPERATURE 300
&THERMOSTAT
TYPE CSVR
&END THERMOSTAT
&END MD
&END MOTION
&MULTIPLE_FORCE_EVALS
FORCE_EVAL_ORDER 2 3
MULTIPLE_SUBSYS F
&END
&FORCE_EVAL
METHOD MIXED
&MIXED
MIXING_TYPE LINEAR_COMBINATION
NGROUPS ${NGROUPS}
&LINEAR
LAMBDA 1.0
MIXED_CDFT TRUE
MIXED_CDFT_COUPLING 1
@IF ( ${USE_DLB} == TRUE )
! Defaults to FALSE
MIXED_CDFT_DLB TRUE
LOAD_SCALE 2
MORE_WORK 0
VERY_OVERLOADED 15
@ENDIF
@IF ( ${COMPUTE_METRIC} == TRUE )
! Defaults to FALSE
MIXED_CDFT_METRIC TRUE
@ENDIF
@IF ( ${WFN_OVERLAP_METHOD} == TRUE )
! Defaults to FALSE
MIXED_CDFT_WFN_OVERLAP TRUE
WFN_RESTART_FILE_NAME HeH-noconstraint-1_0.wfn
@ENDIF
&END LINEAR
&PRINT
&PROGRAM_RUN_INFO
&END
&END PRINT
&END MIXED
@include subsys.inc
&END FORCE_EVAL
&FORCE_EVAL
METHOD QS
@SET BECKE_STR ${BECKE_STR_1}
@SET BECKE_TARGET ${BECKE_TARGET_1}
@SET PROJECT_NAME ${NAME}-state-1
@SET WFN_FILE ${WFN_FILE_1}
@include dft-common-params.inc
&END FORCE_EVAL
&FORCE_EVAL
METHOD QS
@SET BECKE_STR ${BECKE_STR_2}
@SET BECKE_TARGET ${BECKE_TARGET_2}
@SET PROJECT_NAME ${NAME}-state-2
@SET WFN_FILE ${WFN_FILE_2}
@include dft-common-params.inc
&END FORCE_EVAL

6
tests/QS/regtest-cdft-3/HeH-mixed-cdft-2.inp
View file

@ -35,7 +35,9 @@
@SET USE_DLB FALSE
@SET COMPUTE_METRIC TRUE
@SET WFN_OVERLAP_METHOD FALSE
@SET WFN_OVERLAP_METHOD FALSE
@SET NGROUPS 2
&GLOBAL
PROJECT ${PROJECT_NAME}
@ -51,7 +53,7 @@
METHOD MIXED
&MIXED
MIXING_TYPE LINEAR_COMBINATION
NGROUPS 2
NGROUPS ${NGROUPS}
&LINEAR
LAMBDA 1.0
MIXED_CDFT TRUE

4
tests/QS/regtest-cdft-3/HeH-mixed-cdft-3.inp
View file

@ -37,6 +37,8 @@
@SET COMPUTE_METRIC FALSE
@SET WFN_OVERLAP_METHOD TRUE
@SET NGROUPS 2
&GLOBAL
PROJECT ${PROJECT_NAME}
RUN_TYPE ENERGY
@ -51,7 +53,7 @@
METHOD MIXED
&MIXED
MIXING_TYPE LINEAR_COMBINATION
NGROUPS 2
NGROUPS ${NGROUPS}
&LINEAR
LAMBDA 1.0
MIXED_CDFT TRUE

6
tests/QS/regtest-cdft-3/HeH-mixed-cdft-4.inp
View file

@ -35,7 +35,9 @@
@SET USE_DLB FALSE
@SET COMPUTE_METRIC FALSE
@SET WFN_OVERLAP_METHOD FALSE
@SET WFN_OVERLAP_METHOD FALSE
@SET NGROUPS 2
&GLOBAL
PROJECT ${PROJECT_NAME}
@ -54,7 +56,7 @@
METHOD MIXED
&MIXED
MIXING_TYPE LINEAR_COMBINATION
NGROUPS 2
NGROUPS ${NGROUPS}
&LINEAR
LAMBDA 1.0
MIXED_CDFT TRUE

6
tests/QS/regtest-cdft-3/HeH-mixed-cdft-5.inp
View file

@ -35,7 +35,9 @@
@SET USE_DLB TRUE
@SET COMPUTE_METRIC FALSE
@SET WFN_OVERLAP_METHOD FALSE
@SET WFN_OVERLAP_METHOD FALSE
@SET NGROUPS 2
&GLOBAL
PROJECT ${PROJECT_NAME}
@ -54,7 +56,7 @@
METHOD MIXED
&MIXED
MIXING_TYPE LINEAR_COMBINATION
NGROUPS 2
NGROUPS ${NGROUPS}
&LINEAR
LAMBDA 1.0
MIXED_CDFT TRUE

6
tests/QS/regtest-cdft-3/HeH-mixed-cdft-6.inp
View file

@ -35,7 +35,9 @@
@SET USE_DLB FALSE
@SET COMPUTE_METRIC FALSE
@SET WFN_OVERLAP_METHOD FALSE
@SET WFN_OVERLAP_METHOD FALSE
@SET NGROUPS 2
&GLOBAL
PROJECT ${PROJECT_NAME}
@ -65,7 +67,7 @@
METHOD MIXED
&MIXED
MIXING_TYPE LINEAR_COMBINATION
NGROUPS 2
NGROUPS ${NGROUPS}
&LINEAR
LAMBDA 1.0
MIXED_CDFT TRUE

6
tests/QS/regtest-cdft-3/HeH-mixed-cdft-7.inp
View file

@ -35,7 +35,9 @@
@SET USE_DLB TRUE
@SET COMPUTE_METRIC FALSE
@SET WFN_OVERLAP_METHOD FALSE
@SET WFN_OVERLAP_METHOD FALSE
@SET NGROUPS 2
&GLOBAL
PROJECT ${PROJECT_NAME}
@ -65,7 +67,7 @@
METHOD MIXED
&MIXED
MIXING_TYPE LINEAR_COMBINATION
NGROUPS 2
NGROUPS ${NGROUPS}
&LINEAR
LAMBDA 1.0
MIXED_CDFT TRUE

102
tests/QS/regtest-cdft-3/HeH-mixed-cdft-8.inp Normal file
View file

@ -0,0 +1,102 @@
@SET RESTART_WFN TRUE
@SET WFN_FILE_1 HeH-cdft-state-1-1_0.wfn
@SET WFN_FILE_2 HeH-cdft-state-2-1_0.wfn
@SET PROJECT_NAME HeH-mixed-cdft-8
@SET NAME ${PROJECT_NAME}
@SET WRITE_WFN 0
@SET CHARGE 1
@SET WRITE_CUBE FALSE
@SET XYZFILE HeH.xyz
@SET BECKE_ACTIVE TRUE
@SET BECKE_FRAGMENT FALSE
@SET MAX_SCF 5
! He+ H
@SET BECKE_TARGET_1 0.0
@SET BECKE_STR_1 0.707711420058
! He H+
@SET BECKE_TARGET_2 2.0
@SET BECKE_STR_2 -1.008285468010
@SET BECKE_GLOBAL_CUTOFF FALSE
@SET BECKE_CUTOFF_ELEMENT TRUE
@SET BECKE_ADJUST_SIZE TRUE
@SET BECKE_ATOMIC_CHARGES FALSE
@SET BECKE_CAVITY_CONFINE TRUE
@SET BECKE_CAVITY_SHAPE VDW
@SET BECKE_CAVITY_PRINT FALSE
@SET BECKE_SHOULD_SKIP TRUE
@SET BECKE_IN_MEMORY TRUE
@SET USE_DLB FALSE
@SET COMPUTE_METRIC TRUE
@SET WFN_OVERLAP_METHOD TRUE
@SET NGROUPS 1
&GLOBAL
PROJECT ${PROJECT_NAME}
RUN_TYPE ENERGY
PRINT_LEVEL MEDIUM
&END GLOBAL
&MULTIPLE_FORCE_EVALS
FORCE_EVAL_ORDER 2 3
MULTIPLE_SUBSYS F
&END
&FORCE_EVAL
METHOD MIXED
&MIXED
MIXING_TYPE LINEAR_COMBINATION
NGROUPS ${NGROUPS}
&LINEAR
LAMBDA 1.0
MIXED_CDFT TRUE
MIXED_CDFT_COUPLING 1
@IF ( ${USE_DLB} == TRUE )
! Defaults to FALSE
MIXED_CDFT_DLB TRUE
LOAD_SCALE 2
MORE_WORK 0
VERY_OVERLOADED 15
@ENDIF
@IF ( ${COMPUTE_METRIC} == TRUE )
! Defaults to FALSE
MIXED_CDFT_METRIC TRUE
@ENDIF
@IF ( ${WFN_OVERLAP_METHOD} == TRUE )
! Defaults to FALSE
MIXED_CDFT_WFN_OVERLAP TRUE
WFN_RESTART_FILE_NAME HeH-noconstraint-1_0.wfn
@ENDIF
&END LINEAR
&PRINT
&PROGRAM_RUN_INFO
&END
&END PRINT
&END MIXED
@include subsys.inc
&END FORCE_EVAL
&FORCE_EVAL
METHOD QS
@SET BECKE_STR ${BECKE_STR_1}
@SET BECKE_TARGET ${BECKE_TARGET_1}
@SET PROJECT_NAME ${NAME}-state-1
@SET WFN_FILE ${WFN_FILE_1}
@include dft-common-params.inc
&END FORCE_EVAL
&FORCE_EVAL
METHOD QS
@SET BECKE_STR ${BECKE_STR_2}
@SET BECKE_TARGET ${BECKE_TARGET_2}
@SET PROJECT_NAME ${NAME}-state-2
@SET WFN_FILE ${WFN_FILE_2}
@include dft-common-params.inc
&END FORCE_EVAL

105
tests/QS/regtest-cdft-3/HeH-mixed-cdft-9.inp Normal file
View file

@ -0,0 +1,105 @@
@SET RESTART_WFN TRUE
@SET WFN_FILE_1 HeH-cdft-state-1-1_0.wfn
@SET WFN_FILE_2 HeH-cdft-state-2-1_0.wfn
@SET PROJECT_NAME HeH-mixed-cdft-9
@SET NAME ${PROJECT_NAME}
@SET WRITE_WFN 0
@SET CHARGE 1
@SET WRITE_CUBE FALSE
@SET XYZFILE HeH.xyz
@SET BECKE_ACTIVE TRUE
@SET BECKE_FRAGMENT FALSE
@SET MAX_SCF 5
! He+ H
@SET BECKE_TARGET_1 0.0
@SET BECKE_STR_1 0.707711420058
! He H+
@SET BECKE_TARGET_2 2.0
@SET BECKE_STR_2 -1.008285468010
@SET BECKE_GLOBAL_CUTOFF FALSE
@SET BECKE_CUTOFF_ELEMENT TRUE
@SET BECKE_ADJUST_SIZE TRUE
@SET BECKE_ATOMIC_CHARGES FALSE
@SET BECKE_CAVITY_CONFINE TRUE
@SET BECKE_CAVITY_SHAPE VDW
@SET BECKE_CAVITY_PRINT FALSE
@SET BECKE_SHOULD_SKIP TRUE
@SET BECKE_IN_MEMORY TRUE
@SET USE_DLB FALSE
@SET COMPUTE_METRIC FALSE
@SET WFN_OVERLAP_METHOD FALSE
@SET NGROUPS 1
&GLOBAL
PROJECT ${PROJECT_NAME}
RUN_TYPE ENERGY_FORCE
PRINT_LEVEL MEDIUM
&TIMINGS
THRESHOLD 0.001
&END
&END GLOBAL
&MULTIPLE_FORCE_EVALS
FORCE_EVAL_ORDER 2 3
MULTIPLE_SUBSYS F
&END
&FORCE_EVAL
METHOD MIXED
&MIXED
MIXING_TYPE LINEAR_COMBINATION
NGROUPS ${NGROUPS}
&LINEAR
LAMBDA 1.0
MIXED_CDFT TRUE
MIXED_CDFT_COUPLING 1
@IF ( ${USE_DLB} == TRUE )
! Defaults to FALSE
MIXED_CDFT_DLB TRUE
LOAD_SCALE 2
MORE_WORK 0
VERY_OVERLOADED 15
@ENDIF
@IF ( ${COMPUTE_METRIC} == TRUE )
! Defaults to FALSE
MIXED_CDFT_METRIC TRUE
@ENDIF
@IF ( ${WFN_OVERLAP_METHOD} == TRUE )
! Defaults to FALSE
MIXED_CDFT_WFN_OVERLAP TRUE
WFN_RESTART_FILE_NAME HeH-noconstraint-1_0.wfn
@ENDIF
&END LINEAR
&PRINT
&PROGRAM_RUN_INFO
&END
&END PRINT
&END MIXED
@include subsys.inc
&END FORCE_EVAL
&FORCE_EVAL
METHOD QS
@SET BECKE_STR ${BECKE_STR_1}
@SET BECKE_TARGET ${BECKE_TARGET_1}
@SET PROJECT_NAME ${NAME}-state-1
@SET WFN_FILE ${WFN_FILE_1}
@include dft-common-params.inc
&END FORCE_EVAL
&FORCE_EVAL
METHOD QS
@SET BECKE_STR ${BECKE_STR_2}
@SET BECKE_TARGET ${BECKE_TARGET_2}
@SET PROJECT_NAME ${NAME}-state-2
@SET WFN_FILE ${WFN_FILE_2}
@include dft-common-params.inc
&END FORCE_EVAL

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

@ -15,4 +15,11 @@ HeH-mixed-cdft-5.inp 11 1e-13
# MD tests
HeH-mixed-cdft-6.inp 11 2e-11 -2.350318822849201
HeH-mixed-cdft-7.inp 11 2e-11 -2.350318822849201
# Rerun some tests using serial mode (NGROUPS 1)
# Electronic coupling and related
HeH-mixed-cdft-8.inp 73 5e-10 510.592548695828
# Energy/Force
HeH-mixed-cdft-9.inp 11 1e-13 -2.348554384953704
# MD
HeH-mixed-cdft-10.inp 11 2e-13 -2.350318822849201
#EOF

52
tests/QS/regtest-cdft-4/HeH-cdft-broyden-bt1.inp Normal file
View file

@ -0,0 +1,52 @@
@SET RESTART_WFN TRUE
@SET WFN_FILE HeH-noconstraint-1_0.wfn
@SET PROJECT_NAME HeH-cdft-broyden-bt1
@SET WRITE_WFN 0
@SET CHARGE 1
@SET WRITE_CUBE FALSE
@SET XYZFILE HeH.xyz
@SET BECKE_ACTIVE TRUE
@SET BECKE_FRAGMENT FALSE
@SET MAX_SCF 3
! He+ H
@SET BECKE_TARGET_1 0.0
@SET BECKE_STR_1 0.0
! He H+
@SET BECKE_TARGET_2 2.0
@SET BECKE_STR_2 0.0
@SET BECKE_GLOBAL_CUTOFF TRUE
@SET BECKE_CUTOFF_ELEMENT FALSE
@SET BECKE_ADJUST_SIZE TRUE
@SET BECKE_ATOMIC_CHARGES FALSE
@SET BECKE_CAVITY_CONFINE TRUE
@SET BECKE_CAVITY_SHAPE VDW
@SET BECKE_CAVITY_PRINT FALSE
@SET BECKE_SHOULD_SKIP TRUE
@SET BECKE_IN_MEMORY TRUE
@SET OPTIMIZER BROYDEN
@SET JACOBIAN_TYPE FD1
@SET JACOBIAN_FREQ_1 10
@SET JACOBIAN_FREQ_2 1
@SET BROYDEN_TYPE BT1
&GLOBAL
PROJECT ${PROJECT_NAME}
RUN_TYPE ENERGY
PRINT_LEVEL MEDIUM
&END GLOBAL
&FORCE_EVAL
METHOD QS
@SET BECKE_STR ${BECKE_STR_1}
@SET BECKE_TARGET ${BECKE_TARGET_1}
@include dft-common-params.inc
@include subsys.inc
&END FORCE_EVAL

52
tests/QS/regtest-cdft-4/HeH-cdft-broyden-bt1explicit.inp Normal file
View file

@ -0,0 +1,52 @@
@SET RESTART_WFN TRUE
@SET WFN_FILE HeH-noconstraint-1_0.wfn
@SET PROJECT_NAME HeH-cdft-broyden-bt1explicit
@SET WRITE_WFN 0
@SET CHARGE 1
@SET WRITE_CUBE FALSE
@SET XYZFILE HeH.xyz
@SET BECKE_ACTIVE TRUE
@SET BECKE_FRAGMENT FALSE
@SET MAX_SCF 2
! He+ H
@SET BECKE_TARGET_1 0.0
@SET BECKE_STR_1 0.0
! He H+
@SET BECKE_TARGET_2 2.0
@SET BECKE_STR_2 0.0
@SET BECKE_GLOBAL_CUTOFF TRUE
@SET BECKE_CUTOFF_ELEMENT FALSE
@SET BECKE_ADJUST_SIZE TRUE
@SET BECKE_ATOMIC_CHARGES FALSE
@SET BECKE_CAVITY_CONFINE TRUE
@SET BECKE_CAVITY_SHAPE VDW
@SET BECKE_CAVITY_PRINT FALSE
@SET BECKE_SHOULD_SKIP TRUE
@SET BECKE_IN_MEMORY TRUE
@SET OPTIMIZER BROYDEN
@SET JACOBIAN_TYPE FD1
@SET JACOBIAN_FREQ_1 10
@SET JACOBIAN_FREQ_2 1
@SET BROYDEN_TYPE BT1_EXPLICIT
&GLOBAL
PROJECT ${PROJECT_NAME}
RUN_TYPE ENERGY
PRINT_LEVEL MEDIUM
&END GLOBAL
&FORCE_EVAL
METHOD QS
@SET BECKE_STR ${BECKE_STR_1}
@SET BECKE_TARGET ${BECKE_TARGET_1}
@include dft-common-params.inc
@include subsys.inc
&END FORCE_EVAL

52
tests/QS/regtest-cdft-4/HeH-cdft-broyden-bt2.inp Normal file
View file

@ -0,0 +1,52 @@
@SET RESTART_WFN TRUE
@SET WFN_FILE HeH-noconstraint-1_0.wfn
@SET PROJECT_NAME HeH-cdft-broyden-bt2
@SET WRITE_WFN 0
@SET CHARGE 1
@SET WRITE_CUBE FALSE
@SET XYZFILE HeH.xyz
@SET BECKE_ACTIVE TRUE
@SET BECKE_FRAGMENT FALSE
@SET MAX_SCF 3
! He+ H
@SET BECKE_TARGET_1 0.0
@SET BECKE_STR_1 0.0
! He H+
@SET BECKE_TARGET_2 2.0
@SET BECKE_STR_2 0.0
@SET BECKE_GLOBAL_CUTOFF TRUE
@SET BECKE_CUTOFF_ELEMENT FALSE
@SET BECKE_ADJUST_SIZE TRUE
@SET BECKE_ATOMIC_CHARGES FALSE
@SET BECKE_CAVITY_CONFINE TRUE
@SET BECKE_CAVITY_SHAPE VDW
@SET BECKE_CAVITY_PRINT FALSE
@SET BECKE_SHOULD_SKIP TRUE
@SET BECKE_IN_MEMORY TRUE
@SET OPTIMIZER BROYDEN
@SET JACOBIAN_TYPE FD1
@SET JACOBIAN_FREQ_1 10
@SET JACOBIAN_FREQ_2 1
@SET BROYDEN_TYPE BT2
&GLOBAL
PROJECT ${PROJECT_NAME}
RUN_TYPE ENERGY
PRINT_LEVEL MEDIUM
&END GLOBAL
&FORCE_EVAL
METHOD QS
@SET BECKE_STR ${BECKE_STR_1}
@SET BECKE_TARGET ${BECKE_TARGET_1}
@include dft-common-params.inc
@include subsys.inc
&END FORCE_EVAL

64
tests/QS/regtest-cdft-4/HeH-cdft-md-broyden-5-1.inp Normal file
View file

@ -0,0 +1,64 @@
@SET RESTART_WFN TRUE
@SET WFN_FILE HeH-cdft-newton-fd1-1_0.wfn
@SET PROJECT_NAME HeH-cdft-md-broyden-5-1
@SET WRITE_WFN 0
@SET CHARGE 1
@SET WRITE_CUBE FALSE
@SET XYZFILE HeH.xyz
@SET BECKE_ACTIVE TRUE
@SET BECKE_FRAGMENT FALSE
@SET MAX_SCF 5
! He+ H
@SET BECKE_TARGET_1 0.0
@SET BECKE_STR_1 0.731791194315
! He H+
@SET BECKE_TARGET_2 2.0
@SET BECKE_STR_2 0.0
@SET BECKE_GLOBAL_CUTOFF TRUE
@SET BECKE_CUTOFF_ELEMENT FALSE
@SET BECKE_ADJUST_SIZE TRUE
@SET BECKE_ATOMIC_CHARGES FALSE
@SET BECKE_CAVITY_CONFINE TRUE
@SET BECKE_CAVITY_SHAPE VDW
@SET BECKE_CAVITY_PRINT FALSE
@SET BECKE_SHOULD_SKIP TRUE
@SET BECKE_IN_MEMORY TRUE
@SET OPTIMIZER BROYDEN
@SET JACOBIAN_TYPE FD1
@SET JACOBIAN_FREQ_1 5
@SET JACOBIAN_FREQ_2 1
@SET BROYDEN_TYPE BT1
&GLOBAL
PROJECT ${PROJECT_NAME}
RUN_TYPE MD
PRINT_LEVEL MEDIUM
&END GLOBAL
&MOTION
&MD
ENSEMBLE NVT
STEPS 3
TIMESTEP 0.5
TEMPERATURE 300
&THERMOSTAT
TYPE CSVR
&END THERMOSTAT
&END MD
&END MOTION
&FORCE_EVAL
METHOD QS
@SET BECKE_STR ${BECKE_STR_1}
@SET BECKE_TARGET ${BECKE_TARGET_1}
@include dft-common-params.inc
@include subsys.inc
&END FORCE_EVAL

64
tests/QS/regtest-cdft-4/HeH-cdft-md-broyden-5-5.inp Normal file
View file

@ -0,0 +1,64 @@
@SET RESTART_WFN TRUE
@SET WFN_FILE HeH-cdft-newton-fd1-1_0.wfn
@SET PROJECT_NAME HeH-cdft-md-broyden-5-5
@SET WRITE_WFN 0
@SET CHARGE 1
@SET WRITE_CUBE FALSE
@SET XYZFILE HeH.xyz
@SET BECKE_ACTIVE TRUE
@SET BECKE_FRAGMENT FALSE
@SET MAX_SCF 5
! He+ H
@SET BECKE_TARGET_1 0.0
@SET BECKE_STR_1 0.731791194315
! He H+
@SET BECKE_TARGET_2 2.0
@SET BECKE_STR_2 0.0
@SET BECKE_GLOBAL_CUTOFF TRUE
@SET BECKE_CUTOFF_ELEMENT FALSE
@SET BECKE_ADJUST_SIZE TRUE
@SET BECKE_ATOMIC_CHARGES FALSE
@SET BECKE_CAVITY_CONFINE TRUE
@SET BECKE_CAVITY_SHAPE VDW
@SET BECKE_CAVITY_PRINT FALSE
@SET BECKE_SHOULD_SKIP TRUE
@SET BECKE_IN_MEMORY TRUE
@SET OPTIMIZER BROYDEN
@SET JACOBIAN_TYPE FD1
@SET JACOBIAN_FREQ_1 5
@SET JACOBIAN_FREQ_2 5
@SET BROYDEN_TYPE BT1
&GLOBAL
PROJECT ${PROJECT_NAME}
RUN_TYPE MD
PRINT_LEVEL MEDIUM
&END GLOBAL
&MOTION
&MD
ENSEMBLE NVT
STEPS 3
TIMESTEP 0.5
TEMPERATURE 300
&THERMOSTAT
TYPE CSVR
&END THERMOSTAT
&END MD
&END MOTION
&FORCE_EVAL
METHOD QS
@SET BECKE_STR ${BECKE_STR_1}
@SET BECKE_TARGET ${BECKE_TARGET_1}
@include dft-common-params.inc
@include subsys.inc
&END FORCE_EVAL

64
tests/QS/regtest-cdft-4/HeH-cdft-md-newton-1-1.inp Normal file
View file

@ -0,0 +1,64 @@
@SET RESTART_WFN TRUE
@SET WFN_FILE HeH-cdft-newton-fd1-1_0.wfn
@SET PROJECT_NAME HeH-cdft-md-newton-1-1
@SET WRITE_WFN 0
@SET CHARGE 1
@SET WRITE_CUBE FALSE
@SET XYZFILE HeH.xyz
@SET BECKE_ACTIVE TRUE
@SET BECKE_FRAGMENT FALSE
@SET MAX_SCF 5
! He+ H
@SET BECKE_TARGET_1 0.0
@SET BECKE_STR_1 0.731791194315
! He H+
@SET BECKE_TARGET_2 2.0
@SET BECKE_STR_2 0.0
@SET BECKE_GLOBAL_CUTOFF TRUE
@SET BECKE_CUTOFF_ELEMENT FALSE
@SET BECKE_ADJUST_SIZE TRUE
@SET BECKE_ATOMIC_CHARGES FALSE
@SET BECKE_CAVITY_CONFINE TRUE
@SET BECKE_CAVITY_SHAPE VDW
@SET BECKE_CAVITY_PRINT FALSE
@SET BECKE_SHOULD_SKIP TRUE
@SET BECKE_IN_MEMORY TRUE
@SET OPTIMIZER NEWTON
@SET JACOBIAN_TYPE FD1
@SET JACOBIAN_FREQ_1 1
@SET JACOBIAN_FREQ_2 1
@SET BROYDEN_TYPE BT1
&GLOBAL
PROJECT ${PROJECT_NAME}
RUN_TYPE MD
PRINT_LEVEL MEDIUM
&END GLOBAL
&MOTION
&MD
ENSEMBLE NVT
STEPS 3
TIMESTEP 0.5
TEMPERATURE 300
&THERMOSTAT
TYPE CSVR
&END THERMOSTAT
&END MD
&END MOTION
&FORCE_EVAL
METHOD QS
@SET BECKE_STR ${BECKE_STR_1}
@SET BECKE_TARGET ${BECKE_TARGET_1}
@include dft-common-params.inc
@include subsys.inc
&END FORCE_EVAL

52
tests/QS/regtest-cdft-4/HeH-cdft-newton-fd1.inp Normal file
View file

@ -0,0 +1,52 @@
@SET RESTART_WFN TRUE
@SET WFN_FILE HeH-noconstraint-1_0.wfn
@SET PROJECT_NAME HeH-cdft-newton-fd1
@SET WRITE_WFN 1
@SET CHARGE 1
@SET WRITE_CUBE FALSE
@SET XYZFILE HeH.xyz
@SET BECKE_ACTIVE TRUE
@SET BECKE_FRAGMENT FALSE
@SET MAX_SCF 5
! He+ H
@SET BECKE_TARGET_1 0.0
@SET BECKE_STR_1 0.0
! He H+
@SET BECKE_TARGET_2 2.0
@SET BECKE_STR_2 0.0
@SET BECKE_GLOBAL_CUTOFF TRUE
@SET BECKE_CUTOFF_ELEMENT FALSE
@SET BECKE_ADJUST_SIZE TRUE
@SET BECKE_ATOMIC_CHARGES FALSE
@SET BECKE_CAVITY_CONFINE TRUE
@SET BECKE_CAVITY_SHAPE VDW
@SET BECKE_CAVITY_PRINT FALSE
@SET BECKE_SHOULD_SKIP TRUE
@SET BECKE_IN_MEMORY TRUE
@SET OPTIMIZER NEWTON
@SET JACOBIAN_TYPE FD1
@SET JACOBIAN_FREQ_1 1
@SET JACOBIAN_FREQ_2 1
@SET BROYDEN_TYPE BT1
&GLOBAL
PROJECT ${PROJECT_NAME}
RUN_TYPE ENERGY
PRINT_LEVEL MEDIUM
&END GLOBAL
&FORCE_EVAL
METHOD QS
@SET BECKE_STR ${BECKE_STR_1}
@SET BECKE_TARGET ${BECKE_TARGET_1}
@include dft-common-params.inc
@include subsys.inc
&END FORCE_EVAL

52
tests/QS/regtest-cdft-4/HeH-cdft-newton-fd1backward.inp Normal file
View file

@ -0,0 +1,52 @@
@SET RESTART_WFN TRUE
@SET WFN_FILE HeH-noconstraint-1_0.wfn
@SET PROJECT_NAME HeH-cdft-newton-fd1backward
@SET WRITE_WFN 0
@SET CHARGE 1
@SET WRITE_CUBE FALSE
@SET XYZFILE HeH.xyz
@SET BECKE_ACTIVE TRUE
@SET BECKE_FRAGMENT FALSE
@SET MAX_SCF 1
! He+ H
@SET BECKE_TARGET_1 0.0
@SET BECKE_STR_1 0.0
! He H+
@SET BECKE_TARGET_2 2.0
@SET BECKE_STR_2 0.0
@SET BECKE_GLOBAL_CUTOFF TRUE
@SET BECKE_CUTOFF_ELEMENT FALSE
@SET BECKE_ADJUST_SIZE TRUE
@SET BECKE_ATOMIC_CHARGES FALSE
@SET BECKE_CAVITY_CONFINE TRUE
@SET BECKE_CAVITY_SHAPE VDW
@SET BECKE_CAVITY_PRINT FALSE
@SET BECKE_SHOULD_SKIP TRUE
@SET BECKE_IN_MEMORY TRUE
@SET OPTIMIZER NEWTON
@SET JACOBIAN_TYPE FD1_BACKWARD
@SET JACOBIAN_FREQ_1 1
@SET JACOBIAN_FREQ_2 1
@SET BROYDEN_TYPE BT1
&GLOBAL
PROJECT ${PROJECT_NAME}
RUN_TYPE ENERGY
PRINT_LEVEL MEDIUM
&END GLOBAL
&FORCE_EVAL
METHOD QS
@SET BECKE_STR ${BECKE_STR_1}
@SET BECKE_TARGET ${BECKE_TARGET_1}
@include dft-common-params.inc
@include subsys.inc
&END FORCE_EVAL

52
tests/QS/regtest-cdft-4/HeH-cdft-newton-fd1central.inp Normal file
View file

@ -0,0 +1,52 @@
@SET RESTART_WFN TRUE
@SET WFN_FILE HeH-noconstraint-1_0.wfn
@SET PROJECT_NAME HeH-cdft-newton-fd1central
@SET WRITE_WFN 0
@SET CHARGE 1
@SET WRITE_CUBE FALSE
@SET XYZFILE HeH.xyz
@SET BECKE_ACTIVE TRUE
@SET BECKE_FRAGMENT FALSE
@SET MAX_SCF 1
! He+ H
@SET BECKE_TARGET_1 0.0
@SET BECKE_STR_1 0.0
! He H+
@SET BECKE_TARGET_2 2.0
@SET BECKE_STR_2 0.0
@SET BECKE_GLOBAL_CUTOFF TRUE
@SET BECKE_CUTOFF_ELEMENT FALSE
@SET BECKE_ADJUST_SIZE TRUE
@SET BECKE_ATOMIC_CHARGES FALSE
@SET BECKE_CAVITY_CONFINE TRUE
@SET BECKE_CAVITY_SHAPE VDW
@SET BECKE_CAVITY_PRINT FALSE
@SET BECKE_SHOULD_SKIP TRUE
@SET BECKE_IN_MEMORY TRUE
@SET OPTIMIZER NEWTON
@SET JACOBIAN_TYPE FD1_CENTRAL
@SET JACOBIAN_FREQ_1 1
@SET JACOBIAN_FREQ_2 1
@SET BROYDEN_TYPE BT1
&GLOBAL
PROJECT ${PROJECT_NAME}
RUN_TYPE ENERGY
PRINT_LEVEL MEDIUM
&END GLOBAL
&FORCE_EVAL
METHOD QS
@SET BECKE_STR ${BECKE_STR_1}
@SET BECKE_TARGET ${BECKE_TARGET_1}
@include dft-common-params.inc
@include subsys.inc
&END FORCE_EVAL

52
tests/QS/regtest-cdft-4/HeH-cdft-newton-fd2.inp Normal file
View file

@ -0,0 +1,52 @@
@SET RESTART_WFN TRUE
@SET WFN_FILE HeH-noconstraint-1_0.wfn
@SET PROJECT_NAME HeH-cdft-newton-fd2
@SET WRITE_WFN 0
@SET CHARGE 1
@SET WRITE_CUBE FALSE
@SET XYZFILE HeH.xyz
@SET BECKE_ACTIVE TRUE
@SET BECKE_FRAGMENT FALSE
@SET MAX_SCF 1
! He+ H
@SET BECKE_TARGET_1 0.0
@SET BECKE_STR_1 0.0
! He H+
@SET BECKE_TARGET_2 2.0
@SET BECKE_STR_2 0.0
@SET BECKE_GLOBAL_CUTOFF TRUE
@SET BECKE_CUTOFF_ELEMENT FALSE
@SET BECKE_ADJUST_SIZE TRUE
@SET BECKE_ATOMIC_CHARGES FALSE
@SET BECKE_CAVITY_CONFINE TRUE
@SET BECKE_CAVITY_SHAPE VDW
@SET BECKE_CAVITY_PRINT FALSE
@SET BECKE_SHOULD_SKIP TRUE
@SET BECKE_IN_MEMORY TRUE
@SET OPTIMIZER NEWTON
@SET JACOBIAN_TYPE FD2
@SET JACOBIAN_FREQ_1 1
@SET JACOBIAN_FREQ_2 1
@SET BROYDEN_TYPE BT1
&GLOBAL
PROJECT ${PROJECT_NAME}
RUN_TYPE ENERGY
PRINT_LEVEL MEDIUM
&END GLOBAL
&FORCE_EVAL
METHOD QS
@SET BECKE_STR ${BECKE_STR_1}
@SET BECKE_TARGET ${BECKE_TARGET_1}
@include dft-common-params.inc
@include subsys.inc
&END FORCE_EVAL

52
tests/QS/regtest-cdft-4/HeH-cdft-newton-ls-fd2backward.inp Normal file
View file

@ -0,0 +1,52 @@
@SET RESTART_WFN TRUE
@SET WFN_FILE HeH-noconstraint-1_0.wfn
@SET PROJECT_NAME HeH-cdft-newton-ls-fd2backward
@SET WRITE_WFN 0
@SET CHARGE 1
@SET WRITE_CUBE FALSE
@SET XYZFILE HeH.xyz
@SET BECKE_ACTIVE TRUE
@SET BECKE_FRAGMENT FALSE
@SET MAX_SCF 2
! He+ H
@SET BECKE_TARGET_1 0.0
@SET BECKE_STR_1 0.0
! He H+
@SET BECKE_TARGET_2 2.0
@SET BECKE_STR_2 0.0
@SET BECKE_GLOBAL_CUTOFF TRUE
@SET BECKE_CUTOFF_ELEMENT FALSE
@SET BECKE_ADJUST_SIZE TRUE
@SET BECKE_ATOMIC_CHARGES FALSE
@SET BECKE_CAVITY_CONFINE TRUE
@SET BECKE_CAVITY_SHAPE VDW
@SET BECKE_CAVITY_PRINT FALSE
@SET BECKE_SHOULD_SKIP TRUE
@SET BECKE_IN_MEMORY TRUE
@SET OPTIMIZER NEWTON_LS
@SET JACOBIAN_TYPE FD2_BACKWARD
@SET JACOBIAN_FREQ_1 1
@SET JACOBIAN_FREQ_2 1
@SET BROYDEN_TYPE BT1
&GLOBAL
PROJECT ${PROJECT_NAME}
RUN_TYPE ENERGY
PRINT_LEVEL MEDIUM
&END GLOBAL
&FORCE_EVAL
METHOD QS
@SET BECKE_STR ${BECKE_STR_1}
@SET BECKE_TARGET ${BECKE_TARGET_1}
@include dft-common-params.inc
@include subsys.inc
&END FORCE_EVAL

52
tests/QS/regtest-cdft-4/HeH-cdft-secant.inp Normal file
View file

@ -0,0 +1,52 @@
@SET RESTART_WFN TRUE
@SET WFN_FILE HeH-noconstraint-1_0.wfn
@SET PROJECT_NAME HeH-cdft-secant
@SET WRITE_WFN 0
@SET CHARGE 1
@SET WRITE_CUBE FALSE
@SET XYZFILE HeH.xyz
@SET BECKE_ACTIVE TRUE
@SET BECKE_FRAGMENT FALSE
@SET MAX_SCF 2
! He+ H
@SET BECKE_TARGET_1 0.0
@SET BECKE_STR_1 0.0
! He H+
@SET BECKE_TARGET_2 2.0
@SET BECKE_STR_2 0.0
@SET BECKE_GLOBAL_CUTOFF TRUE
@SET BECKE_CUTOFF_ELEMENT FALSE
@SET BECKE_ADJUST_SIZE TRUE
@SET BECKE_ATOMIC_CHARGES FALSE
@SET BECKE_CAVITY_CONFINE TRUE
@SET BECKE_CAVITY_SHAPE VDW
@SET BECKE_CAVITY_PRINT FALSE
@SET BECKE_SHOULD_SKIP TRUE
@SET BECKE_IN_MEMORY TRUE
@SET OPTIMIZER SECANT
@SET JACOBIAN_TYPE FD1
@SET JACOBIAN_FREQ_1 1
@SET JACOBIAN_FREQ_2 1
@SET BROYDEN_TYPE BT1
&GLOBAL
PROJECT ${PROJECT_NAME}
RUN_TYPE ENERGY
PRINT_LEVEL MEDIUM
&END GLOBAL
&FORCE_EVAL
METHOD QS
@SET BECKE_STR ${BECKE_STR_1}
@SET BECKE_TARGET ${BECKE_TARGET_1}
@include dft-common-params.inc
@include subsys.inc
&END FORCE_EVAL

102
tests/QS/regtest-cdft-4/HeH-mixed-cdft-md-newton-1-1.inp Normal file
View file

@ -0,0 +1,102 @@
@SET RESTART_WFN TRUE
@SET WFN_FILE_1 HeH-cdft-newton-fd1-1_0.wfn
@SET WFN_FILE_2 HeH-cdft-newton-fd1-1_0.wfn
@SET PROJECT_NAME HeH-mixed-cdft-md-newton-5-1
@SET NAME ${PROJECT_NAME}
@SET WRITE_WFN 0
@SET CHARGE 1
@SET WRITE_CUBE FALSE
@SET XYZFILE HeH.xyz
@SET BECKE_ACTIVE TRUE
@SET BECKE_FRAGMENT FALSE
@SET MAX_SCF 5
! He+ H
@SET BECKE_TARGET_1 0.0
@SET BECKE_STR_1 0.731791194315
! He H+
@SET BECKE_TARGET_2 2.0
@SET BECKE_STR_2 0.0
@SET BECKE_GLOBAL_CUTOFF TRUE
@SET BECKE_CUTOFF_ELEMENT FALSE
@SET BECKE_ADJUST_SIZE TRUE
@SET BECKE_ATOMIC_CHARGES FALSE
@SET BECKE_CAVITY_CONFINE TRUE
@SET BECKE_CAVITY_SHAPE VDW
@SET BECKE_CAVITY_PRINT FALSE
@SET BECKE_SHOULD_SKIP TRUE
@SET BECKE_IN_MEMORY TRUE
@SET NGROUPS 1
@SET OPTIMIZER NEWTON
@SET JACOBIAN_TYPE FD1
@SET JACOBIAN_FREQ_1 1
@SET JACOBIAN_FREQ_2 1
@SET BROYDEN_TYPE BT1
&GLOBAL
PROJECT ${PROJECT_NAME}
RUN_TYPE MD
PRINT_LEVEL MEDIUM
&END GLOBAL
&MOTION
&MD
ENSEMBLE NVT
STEPS 1
TIMESTEP 0.5
TEMPERATURE 300
&THERMOSTAT
TYPE CSVR
&END THERMOSTAT
&END MD
&END MOTION
&MULTIPLE_FORCE_EVALS
FORCE_EVAL_ORDER 2 3
MULTIPLE_SUBSYS F
&END
&FORCE_EVAL
METHOD MIXED
&MIXED
MIXING_TYPE LINEAR_COMBINATION
NGROUPS ${NGROUPS}
&LINEAR
LAMBDA 1.0
MIXED_CDFT TRUE
MIXED_CDFT_COUPLING 1
&END LINEAR
&PRINT
&PROGRAM_RUN_INFO
&END
&END PRINT
&END MIXED
@include subsys.inc
&END FORCE_EVAL
&FORCE_EVAL
METHOD QS
@SET BECKE_STR ${BECKE_STR_1}
@SET BECKE_TARGET ${BECKE_TARGET_1}
@SET PROJECT_NAME ${NAME}-state-1
@SET WFN_FILE ${WFN_FILE_1}
@include dft-common-params.inc
&END FORCE_EVAL
&FORCE_EVAL
METHOD QS
@SET BECKE_STR ${BECKE_STR_1}
@SET BECKE_TARGET ${BECKE_TARGET_1}
@SET PROJECT_NAME ${NAME}-state-2
@SET WFN_FILE ${WFN_FILE_1}
@include dft-common-params.inc
&END FORCE_EVAL

44
tests/QS/regtest-cdft-4/HeH-noconstraint.inp Normal file
View file

@ -0,0 +1,44 @@
@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 XYZFILE HeH.xyz
@SET BECKE_ACTIVE FALSE
@SET BECKE_FRAGMENT FALSE
@SET MAX_SCF 0
! He+ H
@SET BECKE_TARGET_1 0.0
@SET BECKE_STR_1 0.0
! He H+
@SET BECKE_TARGET_2 2.0
@SET BECKE_STR_2 0.0
@SET BECKE_GLOBAL_CUTOFF TRUE
@SET BECKE_CUTOFF_ELEMENT FALSE
@SET BECKE_ADJUST_SIZE FALSE
@SET BECKE_ATOMIC_CHARGES FALSE
@SET BECKE_CAVITY_CONFINE FALSE
@SET BECKE_CAVITY_SHAPE VDW
@SET BECKE_CAVITY_PRINT FALSE
@SET BECKE_SHOULD_SKIP FALSE
@SET BECKE_IN_MEMORY FALSE
&GLOBAL
PROJECT ${PROJECT_NAME}
RUN_TYPE ENERGY
PRINT_LEVEL MEDIUM
&END GLOBAL
&FORCE_EVAL
METHOD QS
@include dft-common-params.inc
@include subsys.inc
&END FORCE_EVAL

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

@ -0,0 +1,4 @@
2
He 7.500000 7.500000 7.100000
H 7.500000 7.500000 7.900000

25
tests/QS/regtest-cdft-4/TEST_FILES Normal file
View file

@ -0,0 +1,25 @@
# 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 1 1e-13 -2.92909797857414
# These tests use different optimizers for the CDFT constraint
HeH-cdft-secant.inp 1 1e-09 -2.34923603230477
# Newton: test different finite different schemes
HeH-cdft-newton-fd1.inp 11 1e-09 -2.347860618295545
HeH-cdft-newton-fd1backward.inp 1 3e-06 -2.71413183905801
HeH-cdft-newton-fd2.inp 1 4e-06 -2.62570660673001
HeH-cdft-newton-fd1central.inp 1 3e-06 -2.67640494953380
# Newton: use backtracking line search
HeH-cdft-newton-ls-fd2backward.inp 1 2e-06 -2.36109314362094
# Broyden optimizers
HeH-cdft-broyden-bt1.inp 1 1e-09 -2.34812170647367
HeH-cdft-broyden-bt2.inp 1 1e-09 -2.34797348440589
HeH-cdft-broyden-bt1explicit.inp 1 3e-08 -2.34851012154505
# MD test cases
HeH-cdft-md-newton-1-1.inp 11 1e-09 -2.355043088382532
HeH-cdft-md-broyden-5-1.inp 11 1e-09 -2.355042993224512
HeH-cdft-md-broyden-5-5.inp 11 1e-09 -2.355042993523262
# Mixed CDFT: serial mode MD
HeH-mixed-cdft-md-newton-1-1.inp 11 1e-09 -2.348832130105212
#EOF

3
tests/QS/regtest-cdft-4/TEST_FILES_RESET Normal file
View file

@ -0,0 +1,3 @@
#
# add files to be reset here
#

88
tests/QS/regtest-cdft-4/becke_qs.inc Normal file
View file

@ -0,0 +1,88 @@
&QS
METHOD GPW
EPS_DEFAULT 1.0E-12
MAP_CONSISTENT
EXTRAPOLATION ASPC
EXTRAPOLATION_ORDER 3
&CDFT
TYPE_OF_CONSTRAINT BECKE
&OUTER_SCF ON
EPS_SCF 1.0E-2
TYPE BECKE_CONSTRAINT
OPTIMIZER ${OPTIMIZER}
JACOBIAN_TYPE ${JACOBIAN_TYPE}
BROYDEN_TYPE ${BROYDEN_TYPE}
JACOBIAN_STEP 1.0E-4
JACOBIAN_FREQ ${JACOBIAN_FREQ_1} ${JACOBIAN_FREQ_2}
BISECT_TRUST_COUNT 8
EXTRAPOLATION_ORDER 2
MAX_SCF ${MAX_SCF}
STEP_SIZE -1.0
&END
&END CDFT
&BECKE_RESTRAINT
@IF ( ${BECKE_ADJUST_SIZE} == TRUE )
! Defaults to false
ADJUST_SIZE TRUE
ATOMIC_RADII 0.460 0.320
@ENDIF
@IF ( ${BECKE_ATOMIC_CHARGES} == TRUE )
! Defaults to false
ATOMIC_CHARGES TRUE
@ENDIF
STRENGTH ${BECKE_STR}
! The constraint target: sum_i coeff_i * N_i
! where N_i is the number of VALENCE electrons on i
TARGET ${BECKE_TARGET}
@IF ( ${BECKE_CUTOFF_ELEMENT} == TRUE )
CUTOFF_TYPE ELEMENT
ELEMENT_CUTOFF 2.0 2.0
@ENDIF
@IF ( ${BECKE_GLOBAL_CUTOFF} == TRUE )
CUTOFF_TYPE GLOBAL
GLOBAL_CUTOFF 2.0
@ENDIF
@IF ( ${BECKE_IN_MEMORY} == TRUE )
! Defaults to false
IN_MEMORY TRUE
@ENDIF
@IF ( ${BECKE_CAVITY_CONFINE} == TRUE )
! Defaults to NONE
CONFINE CAVITY
EPS_CAVITY 1.0E-6
CAVITY_SHAPE ${BECKE_CAVITY_SHAPE}
! For shape single
CAVITY_RADIUS 1.3
CAVITY_USE_BOHR FALSE
@ENDIF
@IF ( ${BECKE_SHOULD_SKIP} == TRUE )
! Defaults to false
SHOULD_SKIP TRUE
@ENDIF
@IF ( ${BECKE_CAVITY_PRINT} == TRUE )
! Defaults to false
CAVITY_PRINT TRUE
@ENDIF
@IF ( ${BECKE_FRAGMENT} == TRUE )
! Need absolute weight to use fragment densities
ATOMS 1
COEFF 1
FRAGMENT_DENSITIES
FRAGMENT_A_FILE_NAME He+-noconstraint-ELECTRON_DENSITY-1_0.cube
FRAGMENT_B_FILE_NAME H-noconstraint-ELECTRON_DENSITY-1_0.cube
@ENDIF
@IF ( ${BECKE_FRAGMENT} == FALSE )
ATOMS 1..2
COEFF 1 -1
@ENDIF
CONSTRAINT_TYPE TOTAL
&PROGRAM_RUN_INFO ON
&EACH
QS_SCF 1
&END EACH
COMMON_ITERATION_LEVELS 2
ADD_LAST NUMERIC
FILENAME ./${PROJECT_NAME}
&END PROGRAM_RUN_INFO
&END BECKE_RESTRAINT
&END QS

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

@ -0,0 +1,82 @@
&DFT
@IF ( ${BECKE_ACTIVE} == TRUE )
@include becke_qs.inc
@ENDIF
@IF ( ${BECKE_ACTIVE} == FALSE )
&QS
METHOD GPW
EPS_DEFAULT 1.0E-12
MAP_CONSISTENT
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
&END
&END

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

@ -0,0 +1,20 @@
&SUBSYS
&CELL
PERIODIC XYZ
ABC 15. 15. 15.
&END CELL
&TOPOLOGY
COORD_FILE_FORMAT XYZ
COORD_FILE_NAME ${XYZFILE}
&CENTER_COORDINATES OFF
&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

2
tests/TEST_DIRS
View file

@ -3,6 +3,8 @@
# Directories have been reordered according the execution time needed for a gfortran pdbg run using 2 MPI tasks
# in case a new directory is added just add it at the top of the list..
# the order will be regularly checked and modified...
QMMM/QS/regtest-cdft parallel mpiranks>1
QS/regtest-cdft-4
QS/regtest-cdft-3 parallel mpiranks>1
QS/regtest-cdft-2
QS/regtest-cdft-1

3
tests/TEST_TYPES
View file

@ -1,4 +1,4 @@
74
75
Total energy:!3
POTENTIAL ENERGY!4
Total energy \[eV\]:!4
@ -73,6 +73,7 @@ Current value of constraint !6
SUM OF ATOMIC FORCES !8
Diabatic electronic coupling !5
Diabatic coupling !7
Charge transfer energy !7
#
# these are the tests the can be selected for regtesting.
# do regtest will grep for test_grep (first column) and look if the numeric value