diff --git a/src/force_env_methods.F b/src/force_env_methods.F index 8744a7204a..58c145b113 100644 --- a/src/force_env_methods.F +++ b/src/force_env_methods.F @@ -1743,7 +1743,8 @@ CONTAINS ! Print embedding potential for restart CALL print_embed_restart(force_env%sub_force_env(ref_subsys_number)%force_env%qs_env, & - opt_embed%dimen_aux, opt_embed%embed_pot_coef, embed_pot, i_iter) + opt_embed%dimen_aux, opt_embed%embed_pot_coef, embed_pot, i_iter, & + spin_embed_pot, opt_embed%open_shell_embed) ! Print information and check convergence CALL print_emb_opt_info(output_unit, i_iter, opt_embed) CALL conv_check_embed(opt_embed, diff_rho_r, diff_rho_spin, output_unit) diff --git a/src/optimize_embedding_potential.F b/src/optimize_embedding_potential.F index e3fc62e818..e809c6087a 100644 --- a/src/optimize_embedding_potential.F +++ b/src/optimize_embedding_potential.F @@ -442,7 +442,8 @@ CONTAINS !> \param embed_pot_coef ... !> \param open_shell_embed ... ! ************************************************************************************************** - SUBROUTINE read_embed_pot(qs_env, embed_pot, spin_embed_pot, section, embed_pot_coef, open_shell_embed) + SUBROUTINE read_embed_pot(qs_env, embed_pot, spin_embed_pot, section, embed_pot_coef, & + open_shell_embed) TYPE(qs_environment_type), POINTER :: qs_env TYPE(pw_p_type), POINTER :: embed_pot, spin_embed_pot TYPE(section_vals_type), POINTER :: section @@ -451,7 +452,8 @@ CONTAINS CHARACTER(LEN=default_path_length) :: filename INTEGER :: dimen_aux, dimen_restart_basis, & - l_global, LLL, nrow_local, restart_unit + dimen_var_aux, l_global, LLL, & + nrow_local, restart_unit INTEGER, DIMENSION(:), POINTER :: row_indices REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: coef, coef_read TYPE(cp_blacs_env_type), POINTER :: blacs_env @@ -461,6 +463,11 @@ CONTAINS ! Get the vector dimension CALL find_aux_dimen(qs_env, dimen_aux) + IF (open_shell_embed) THEN + dimen_var_aux = dimen_aux*2 + ELSE + dimen_var_aux = dimen_aux + ENDIF ! We need a temporary vector of coefficients CALL get_qs_env(qs_env=qs_env, & @@ -469,13 +476,13 @@ CONTAINS NULLIFY (my_embed_pot_coef, fm_struct) CALL cp_blacs_env_create(blacs_env=blacs_env, para_env=para_env) CALL cp_fm_struct_create(fm_struct, para_env=para_env, context=blacs_env, & - nrow_global=dimen_aux, ncol_global=1) + nrow_global=dimen_var_aux, ncol_global=1) CALL cp_fm_create(my_embed_pot_coef, fm_struct, name="my_pot_coef") IF (.NOT. (PRESENT(embed_pot_coef))) THEN NULLIFY (embed_pot_coef) CALL cp_fm_create(my_embed_pot_coef, fm_struct, name="pot_coef") - CALL cp_fm_set_all(my_embed_pot_coef, 0.0_dp) ENDIF + CALL cp_fm_struct_release(fm_struct) CALL cp_fm_set_all(my_embed_pot_coef, 0.0_dp) @@ -483,13 +490,11 @@ CONTAINS restart_unit = -1 ! Allocate the attay to read the coefficients - ALLOCATE (coef(dimen_aux)) + ALLOCATE (coef(dimen_var_aux)) coef = 0.0_dp IF (para_env%ionode) THEN - ALLOCATE (coef_read(dimen_aux)) - coef_read = 0.0_dp ! Get the restart file name CALL embed_restart_file_name(filename, section) @@ -504,6 +509,9 @@ CONTAINS IF (.NOT. (dimen_restart_basis == dimen_aux)) & CPABORT("Wrong dimension of the embedding basis in the restart file.") + ALLOCATE (coef_read(dimen_var_aux)) + coef_read = 0.0_dp + READ (restart_unit) coef_read coef(:) = coef_read(:) DEALLOCATE (coef_read) @@ -526,6 +534,7 @@ CONTAINS l_global = row_indices(LLL) my_embed_pot_coef%local_data(LLL, 1) = coef(l_global) ENDDO + DEALLOCATE (coef) ! Copy to the my_embed_pot_coef to embed_pot_coef @@ -1894,13 +1903,18 @@ CONTAINS !> \param embed_pot_coef ... !> \param embed_pot ... !> \param i_iter ... +!> \param embed_pot_spin ... +!> \param open_shell_embed ... ! ************************************************************************************************** - SUBROUTINE print_embed_restart(qs_env, dimen_aux, embed_pot_coef, embed_pot, i_iter) + SUBROUTINE print_embed_restart(qs_env, dimen_aux, embed_pot_coef, embed_pot, i_iter, & + embed_pot_spin, open_shell_embed) TYPE(qs_environment_type), POINTER :: qs_env INTEGER :: dimen_aux TYPE(cp_fm_type), POINTER :: embed_pot_coef TYPE(pw_p_type), POINTER :: embed_pot INTEGER :: i_iter + TYPE(pw_p_type), POINTER :: embed_pot_spin + LOGICAL :: open_shell_embed CHARACTER(LEN=default_path_length) :: filename, my_pos_cube, title INTEGER :: unit_nr @@ -1930,7 +1944,7 @@ CONTAINS ENDIF ENDIF - ! Second a cube file + ! Second, cube files dft_section => section_vals_get_subs_vals(input, "DFT") CALL qs_subsys_get(subsys, particles=particles) @@ -1948,6 +1962,19 @@ CONTAINS stride=section_get_ivals(dft_section, "QS%OPT_EMBED%EMBED_POT_CUBE%STRIDE")) CALL cp_print_key_finished_output(unit_nr, logger, input, & "DFT%QS%OPT_EMBED%EMBED_POT_CUBE") + IF (open_shell_embed) THEN ! Print spin part of the embedding potential + my_pos_cube = "REWIND" + WRITE (filename, '(a15,I3.3)') "spin_embed_pot_", i_iter + unit_nr = cp_print_key_unit_nr(logger, input, "DFT%QS%OPT_EMBED%EMBED_POT_CUBE", & + extension=".cube", middle_name=TRIM(filename), file_position=my_pos_cube, & + log_filename=.FALSE.) + + WRITE (title, *) "SPIN EMBEDDING POTENTIAL at optimization step ", i_iter + CALL cp_pw_to_cube(embed_pot_spin%pw, unit_nr, title, particles=particles, & + stride=section_get_ivals(dft_section, "QS%OPT_EMBED%EMBED_POT_CUBE%STRIDE")) + CALL cp_print_key_finished_output(unit_nr, logger, input, & + "DFT%QS%OPT_EMBED%EMBED_POT_CUBE") + ENDIF ENDIF END SUBROUTINE print_embed_restart diff --git a/tests/QS/regtest-embed/H_H_pbe_pbe0_triplet_restart.inp b/tests/QS/regtest-embed/H_H_pbe_pbe0_triplet_restart.inp new file mode 100644 index 0000000000..9fa749f436 --- /dev/null +++ b/tests/QS/regtest-embed/H_H_pbe_pbe0_triplet_restart.inp @@ -0,0 +1,312 @@ +#CPQA DEPENDS H_H_pbe_pbe0_triplet.inp +! +! Test restart in the spin-unrestricted case +! +&GLOBAL + PROJECT h_h_pbe_pbe0_triplet_restart + PRINT_LEVEL MEDIUM + RUN_TYPE ENERGY +&END GLOBAL +&MULTIPLE_FORCE_EVALS + FORCE_EVAL_ORDER 2 3 4 5 + MULTIPLE_SUBSYS T +&END +&FORCE_EVAL + METHOD EMBED + &EMBED + NGROUPS 1 + &MAPPING + &FORCE_EVAL_EMBED + &FRAGMENT 1 + 1 1 + &END + &FRAGMENT 2 + 2 2 + &END + &FRAGMENT 3 + 1 2 + &END + &END + &FORCE_EVAL 1 + &FRAGMENT 1 + 1 1 + MAP 1 + &END + &END + &FORCE_EVAL 2 + &FRAGMENT 1 + 1 1 + MAP 2 + &END + &END + &FORCE_EVAL 3 + &FRAGMENT 1 + 1 2 + MAP 3 + &END + &END + &FORCE_EVAL 4 + &FRAGMENT 1 + 1 1 + MAP 2 + &END + &END + &END + &END EMBED + &SUBSYS + &CELL + ABC [angstrom] 5.000 5.000 5.000 + &END CELL + &KIND H + BASIS_SET cc-TZ + RI_AUX_BASIS_SET RI_TZ + POTENTIAL GTH-HF-q1 + &END KIND + &KIND O + BASIS_SET cc-TZ + RI_AUX_BASIS_SET RI_TZ + POTENTIAL GTH-HF-q6 + &END KIND + &COORD + H 1.75 2.75 0.0 + H 1.75 4.25 0.0 + &END + &END SUBSYS +&END FORCE_EVAL + +! Subsys 1 + +&FORCE_EVAL + METHOD Quickstep + &DFT + BASIS_SET_FILE_NAME BASIS_RI_cc-TZ + POTENTIAL_FILE_NAME HF_POTENTIALS + UKS .TRUE. + MULTIPLICITY 2 + &MGRID + CUTOFF 100 + REL_CUTOFF 20 + &END MGRID + &POISSON + &END POISSON + &QS + METHOD GPW + EPS_DEFAULT 1.0E-15 + EPS_PGF_ORB 1.0E-30 + &END QS + &SCF + &OT + PRECONDITIONER FULL_ALL + &END + SCF_GUESS ATOMIC + MAX_SCF 100 + &PRINT + &RESTART OFF + &END + &END + &END SCF + &XC + &XC_FUNCTIONAL PBE + &END XC_FUNCTIONAL + &END XC + &END DFT + &SUBSYS + &CELL + ABC [angstrom] 5.000 5.000 5.000 + &END CELL + &KIND H + BASIS_SET cc-TZ + RI_AUX_BASIS_SET RI_TZ + POTENTIAL GTH-HF-q1 + &END KIND + &COORD + H 1.75 2.75 0.0 + &END + &END SUBSYS +&END FORCE_EVAL + +! Subsys 2 + +&FORCE_EVAL + METHOD Quickstep + &DFT + BASIS_SET_FILE_NAME BASIS_RI_cc-TZ + POTENTIAL_FILE_NAME HF_POTENTIALS + UKS .TRUE. + MULTIPLICITY 2 + &MGRID + CUTOFF 100 + REL_CUTOFF 20 + &END MGRID + &POISSON + &END POISSON + &QS + CLUSTER_EMBED_SUBSYS .TRUE. + METHOD GPW + EPS_DEFAULT 1.0E-15 + EPS_PGF_ORB 1.0E-30 + &END QS + &SCF + &OT + PRECONDITIONER FULL_ALL + &END + SCF_GUESS ATOMIC + MAX_SCF 100 + &PRINT + &RESTART OFF + &END + &END + &END SCF + &XC + &XC_FUNCTIONAL PBE + &END XC_FUNCTIONAL + &END XC + &END DFT + &SUBSYS + &CELL + ABC [angstrom] 5.000 5.000 5.000 + &END CELL + &KIND H + BASIS_SET cc-TZ + RI_AUX_BASIS_SET RI_TZ + POTENTIAL GTH-HF-q1 + &END KIND + &COORD + H 1.75 4.25 0.0 + &END + &END SUBSYS +&END FORCE_EVAL + +! Total system + +&FORCE_EVAL + METHOD Quickstep + &DFT + UKS .TRUE. + MULTIPLICITY 3 + BASIS_SET_FILE_NAME BASIS_RI_cc-TZ + POTENTIAL_FILE_NAME HF_POTENTIALS + &PRINT + &E_DENSITY_CUBE HIGH + &END + &END + &MGRID + CUTOFF 100 + REL_CUTOFF 20 + &END MGRID + &POISSON + &END POISSON + &QS + REF_EMBED_SUBSYS .TRUE. + METHOD GPW + EPS_DEFAULT 1.0E-15 + EPS_PGF_ORB 1.0E-30 + &OPT_EMBED + REG_LAMBDA 0.00001 + N_ITER 50 + DENS_CONV_INT 0.5 + SPIN_DENS_CONV_INT 0.5 + DENS_CONV_MAX 0.025 + &EMBED_POT_CUBE MEDIUM + &END + READ_EMBED_POT .TRUE. + EMBED_RESTART_FILE_NAME h_h_pbe_pbe0_triplet-embed_pot_002-1_0.wfn + &END + &END QS + &SCF + &OT + PRECONDITIONER FULL_ALL + &END + SCF_GUESS ATOMIC + MAX_SCF 100 + &PRINT + &RESTART OFF + &END + &END + &END SCF + &XC + &XC_FUNCTIONAL PBE + &END XC_FUNCTIONAL + &END XC + &END DFT + &SUBSYS + &CELL + ABC [angstrom] 5.000 5.000 5.000 + &END CELL + &KIND H + BASIS_SET cc-TZ + RI_AUX_BASIS_SET RI_TZ + POTENTIAL GTH-HF-q1 + &END KIND + &COORD + H 1.75 2.75 0.0 + H 1.75 4.25 0.0 + &END + &END SUBSYS +&END FORCE_EVAL + +! Higher level calculation on subsys 2 + +&FORCE_EVAL + METHOD Quickstep + &DFT + UKS .TRUE. + MULTIPLICITY 2 + BASIS_SET_FILE_NAME BASIS_RI_cc-TZ + POTENTIAL_FILE_NAME HF_POTENTIALS + &MGRID + CUTOFF 100 + REL_CUTOFF 20 + &END MGRID + &POISSON + &END POISSON + &QS + HIGH_LEVEL_EMBED_SUBSYS + METHOD GPW + EPS_DEFAULT 1.0E-15 + EPS_PGF_ORB 1.0E-30 + &END QS + &SCF + &OT + PRECONDITIONER FULL_ALL + &END + SCF_GUESS ATOMIC + MAX_SCF 100 + &PRINT + &RESTART OFF + &END + &END + &END SCF + &XC + &XC_FUNCTIONAL PBE + &PBE + SCALE_X 0.75 + SCALE_C 1.0 + &END + &END XC_FUNCTIONAL + &HF + FRACTION 0.25 + &INTERACTION_POTENTIAL + POTENTIAL_TYPE TRUNCATED + CUTOFF_RADIUS 2.45 + T_C_G_DATA t_c_g.dat + &END + + &END + &END XC + &END DFT + &SUBSYS + &CELL + ABC [angstrom] 5.000 5.000 5.000 + &END CELL + &KIND H + BASIS_SET cc-TZ + RI_AUX_BASIS_SET RI_TZ + POTENTIAL GTH-HF-q1 + &END KIND + &COORD + H 1.75 4.25 0.0 + &END + &END SUBSYS +&END FORCE_EVAL + diff --git a/tests/QS/regtest-embed/TEST_FILES b/tests/QS/regtest-embed/TEST_FILES index 6b7edb1069..e6425d6a02 100644 --- a/tests/QS/regtest-embed/TEST_FILES +++ b/tests/QS/regtest-embed/TEST_FILES @@ -3,5 +3,6 @@ H2O_H2_pbe_mp2.inp 11 1e-9 -1 H2O_H2_pbe_rpa_restart.inp 11 2e-8 -18.3610866706 H4_H8_pbe_pbe0_const_pot.inp 11 1e-9 -3.97327602819 H_H_pbe_pbe0_triplet.inp 11 1e-9 -0.94029854921 +H_H_pbe_pbe0_triplet_restart.inp 11 1e-9 -0.94029854921 H_H_pbe_pbe0_singlet_roks.inp 11 1e-9 -1.05060718524 #EOF