diff --git a/src/bsse.F b/src/bsse.F index 54243ecded..3abfc8c654 100644 --- a/src/bsse.F +++ b/src/bsse.F @@ -13,7 +13,6 @@ ! ************************************************************************************************** MODULE bsse USE atomic_kind_types, ONLY: get_atomic_kind - USE cell_types, ONLY: cell_type USE cp2k_info, ONLY: write_restart_header USE cp_external_control, ONLY: external_control USE cp_log_handling, ONLY: cp_get_default_logger,& @@ -236,7 +235,6 @@ CONTAINS my_targ, n_rep, num_of_frag, old_size, & present_charge, present_multpl INTEGER, DIMENSION(:), POINTER :: atom_index, atom_list, my_conf, tmplist - TYPE(cell_type), POINTER :: cell TYPE(cp_subsys_type), POINTER :: subsys TYPE(mp_para_env_type), POINTER :: para_env TYPE(particle_list_type), POINTER :: particles @@ -245,7 +243,7 @@ CONTAINS CALL section_vals_get(n_frags, n_repetition=num_of_frag) CPASSERT(SIZE(conf) == num_of_frag) - NULLIFY (subsys, particles, para_env, cell, atom_index, atom_type, tmplist, & + NULLIFY (subsys, particles, para_env, atom_index, atom_type, tmplist, & force_env_section) CALL force_env_get(force_env, force_env_section=force_env_section) CALL section_vals_val_get(force_env_section, "METHOD", i_val=method_name_id) @@ -255,8 +253,7 @@ CONTAINS ALLOCATE (my_conf(SIZE(conf))) my_conf = conf - CALL force_env_get(force_env=force_env, subsys=subsys, para_env=para_env, & - cell=cell) + CALL force_env_get(force_env=force_env, subsys=subsys, para_env=para_env) CALL cp_subsys_get(subsys, particles=particles) isize = 0 ALLOCATE (atom_index(isize)) @@ -313,7 +310,7 @@ CONTAINS TYPE(cp_subsys_type), POINTER :: subsys_loc NULLIFY (subsys_loc) CALL create_small_subsys(subsys_loc, big_subsys=subsys, & - small_para_env=para_env, small_cell=cell, sub_atom_index=atom_index, & + small_para_env=para_env, small_cell=subsys%cell, sub_atom_index=atom_index, & sub_atom_kind_name=atom_type, para_env=para_env, & force_env_section=force_env_section, subsys_section=subsys_section) @@ -334,7 +331,7 @@ CONTAINS DEALLOCATE (qs_env) END BLOCK ELSE - CPABORT("") + CPABORT("BSSE calculation is only available for QuickStep method!") END IF DEALLOCATE (atom_index) DEALLOCATE (atom_type) diff --git a/src/cp_ddapc_types.F b/src/cp_ddapc_types.F index eb1f282abd..9ae73b0dc7 100644 --- a/src/cp_ddapc_types.F +++ b/src/cp_ddapc_types.F @@ -18,7 +18,6 @@ MODULE cp_ddapc_types ewald_ddapc_pot,& solvation_ddapc_pot USE cp_log_handling, ONLY: cp_get_default_logger,& - cp_logger_get_default_io_unit,& cp_logger_type USE cp_output_handling, ONLY: cp_printkey_is_on USE ewald_spline_util, ONLY: Setup_Ewald_Spline @@ -28,6 +27,7 @@ MODULE cp_ddapc_types section_vals_val_get USE kinds, ONLY: dp USE mathconstants, ONLY: pi + USE mathlib, ONLY: det_3x3 USE message_passing, ONLY: mp_para_env_type USE particle_types, ONLY: particle_type USE pw_grid_types, ONLY: pw_grid_type @@ -224,7 +224,7 @@ CONTAINS TYPE(section_vals_type), POINTER :: force_env_section, subsys_section TYPE(mp_para_env_type), POINTER :: para_env - INTEGER :: my_val, npts(3), output_unit + INTEGER :: my_val, npts(3) INTEGER, DIMENSION(:), POINTER :: ngrids LOGICAL :: analyt, decoupling, & do_qmmm_periodic_decpl, do_restraint, & @@ -238,7 +238,6 @@ CONTAINS solvation_section logger => cp_get_default_logger() - output_unit = cp_logger_get_default_io_unit(logger) CPASSERT(.NOT. ASSOCIATED(cp_ddapc_ewald)) ALLOCATE (cp_ddapc_ewald) NULLIFY (cp_ddapc_ewald%pw_grid_mm, & @@ -268,12 +267,10 @@ CONTAINS cp_ddapc_ewald%do_restraint = do_restraint .OR. do_restraintB ! Determining the tasks and further check IF (do_qmmm_periodic_decpl .AND. decoupling) THEN - ! Check than an additional POISSON section has not been defined. In case write a warning - IF (output_unit > 0) & - WRITE (output_unit, '(T2,"WARNING",A)') & - "A calculation with the QMMM periodic model has been requested.", & - "The explicit POISSON section in DFT section will be IGNORED.", & - "QM Electrostatic controlled only by the PERIODIC section in QMMM section" + CALL cp_warn(__LOCATION__, & + "A calculation with the QMMM periodic model has been requested. "// & + "The explicit POISSON section in DFT section will be IGNORED. "// & + "QM Electrostatic controlled only by the PERIODIC section in QMMM section") decoupling = .FALSE. END IF IF (decoupling) THEN @@ -301,7 +298,7 @@ CONTAINS NULLIFY (LG, gx, gy, gz) hmat = qm_cell%hmat - CALL eval_lg(multipole_section, hmat, qm_cell%deth, LG, gx, gy, gz) + CALL eval_lg(multipole_section, hmat, LG, gx, gy, gz) grid_print_section => section_vals_get_subs_vals(force_env_section, "PRINT%GRID_INFORMATION") CALL Setup_Ewald_Spline(pw_grid=cp_ddapc_ewald%pw_grid_qm, pw_pool=cp_ddapc_ewald%pw_pool_qm, & coeff=cp_ddapc_ewald%coeff_qm, LG=LG, gx=gx, gy=gy, gz=gz, hmat=hmat, npts=npts, & @@ -316,7 +313,7 @@ CONTAINS cell_section => section_vals_get_subs_vals(subsys_section, "CELL") CALL read_cell(mm_cell, dummy_cell, cell_section=cell_section, para_env=para_env) hmat = mm_cell%hmat - CALL eval_lg(multipole_section, hmat, mm_cell%deth, LG, gx, gy, gz) + CALL eval_lg(multipole_section, hmat, LG, gx, gy, gz) grid_print_section => section_vals_get_subs_vals(force_env_section, "PRINT%GRID_INFORMATION") CALL Setup_Ewald_Spline(pw_grid=cp_ddapc_ewald%pw_grid_mm, pw_pool=cp_ddapc_ewald%pw_pool_mm, & coeff=cp_ddapc_ewald%coeff_mm, LG=LG, gx=gx, gy=gy, gz=gz, hmat=hmat, npts=npts, & @@ -337,7 +334,6 @@ CONTAINS !> \brief ... !> \param multipole_section ... !> \param hmat ... -!> \param deth ... !> \param LG ... !> \param gx ... !> \param gy ... @@ -346,16 +342,17 @@ CONTAINS !> none !> \author Teodoro Laino - [tlaino] ! ************************************************************************************************** - SUBROUTINE eval_lg(multipole_section, hmat, deth, LG, gx, gy, gz) + SUBROUTINE eval_lg(multipole_section, hmat, LG, gx, gy, gz) TYPE(section_vals_type), POINTER :: multipole_section - REAL(KIND=dp), INTENT(IN) :: hmat(3, 3), deth + REAL(KIND=dp), INTENT(IN) :: hmat(3, 3) REAL(KIND=dp), DIMENSION(:), POINTER :: LG, gx, gy, gz INTEGER :: i, k1, k2, k3, n_rep, ndim, nmax1, & nmax2, nmax3 - REAL(KIND=dp) :: alpha, eps, fac, fs, fvec(3), galpha, & - gsq, gsqi, rcut, tol, tol1 + REAL(KIND=dp) :: alpha, deth, eps, fac, fs, fvec(3), & + galpha, gsq, gsqi, rcut, tol, tol1 + deth = ABS(det_3x3(hmat)) rcut = MIN(hmat(1, 1), hmat(2, 2), hmat(3, 3))/2.0_dp CALL section_vals_val_get(multipole_section, "RCUT", n_rep_val=n_rep) IF (n_rep == 1) CALL section_vals_val_get(multipole_section, "RCUT", r_val=rcut) diff --git a/src/cp_subsys_methods.F b/src/cp_subsys_methods.F index 9d5030b1e0..12bea09964 100644 --- a/src/cp_subsys_methods.F +++ b/src/cp_subsys_methods.F @@ -18,6 +18,7 @@ MODULE cp_subsys_methods atomic_kind_list_type USE atomic_kind_types, ONLY: atomic_kind_type USE atprop_types, ONLY: atprop_create + USE cell_methods, ONLY: write_cell USE cell_types, ONLY: cell_retain,& cell_type USE colvar_methods, ONLY: colvar_read @@ -237,7 +238,7 @@ CONTAINS !> \brief updates the molecule information of the given subsys !> \param small_subsys the subsys to create !> \param big_subsys the superset of small_subsys -!> \param small_cell ... +!> \param small_cell the cell of small_subsys !> \param small_para_env the parallel environment for the new (small) !> subsys !> \param sub_atom_index indexes of the atoms that should be in small_subsys @@ -252,6 +253,9 @@ CONTAINS !> \note !> not really ready to be used with different para_envs for the small !> and big part +!> qmmm_env_create() should be the only instance where this subroutine is called +!> with small_cell distinct from big_subsys%cell; otherwise it could be possible +!> to drop the small_cell as input parameter entirely. ! ************************************************************************************************** SUBROUTINE create_small_subsys(small_subsys, big_subsys, small_cell, & small_para_env, sub_atom_index, sub_atom_kind_name, & @@ -299,6 +303,8 @@ CONTAINS !----------------------------------------------------------------------------- topology%cell => small_cell CALL cell_retain(small_cell) + CALL write_cell(small_cell, subsys_section, tag="CELL") + CALL write_cell(small_cell, subsys_section, tag="CELL_REF") !----------------------------------------------------------------------------- !----------------------------------------------------------------------------- @@ -378,7 +384,8 @@ CONTAINS CALL molecule_list_create(mols, els_ptr=molecule_set) CALL molecule_kind_list_create(mol_kinds, els_ptr=molecule_kind_set) CALL cp_subsys_set(small_subsys, particles=particles, atomic_kinds=atomic_kinds, & - molecules=mols, molecule_kinds=mol_kinds) + molecules=mols, molecule_kinds=mol_kinds, cell=small_cell, & + cell_ref=small_cell, use_ref_cell=.FALSE.) CALL particle_list_release(particles) CALL atomic_kind_list_release(atomic_kinds) CALL molecule_list_release(mols) diff --git a/src/input_cp2k_subsys.F b/src/input_cp2k_subsys.F index 1694fe9168..f57d48f45b 100644 --- a/src/input_cp2k_subsys.F +++ b/src/input_cp2k_subsys.F @@ -87,15 +87,15 @@ CONTAINS "2. The lengths and angles of cell vectors set by ABC and ALPHA_BETA_GAMMA;"//newline// & "3. The vectors set by A, B, C together;"//newline// & "4. If none above exist, the external file set by TOPOLOGY/COORD_FILE_NAME with "// & - "certain TOPOLOGY/COORD_FILE_FORMAT may also be parsed. This behavior has NOT "// & - "been fully tested for BSSE as well as QM/MM calculations.") + "suitable TOPOLOGY/COORD_FILE_FORMAT may also be parsed for FORCE_EVAL/SUBSYS/CELL "// & + "but not for FORCE_EVAL/QMMM/CELL.") CALL create_cell_section_low(section, periodic) NULLIFY (subsection) CALL section_create(subsection, __LOCATION__, "CELL_REF", & - description="Input parameters needed to set up the reference cell. "// & - "This option can be used to keep the FFT grid fixed while "// & - "running a cell optimization or NpT molecular dynamics. "// & + description="Input parameters needed to set up the reference cell for "// & + "FORCE_EVAL/SUBSYS/CELL. This option can be used to keep the FFT grid "// & + "fixed while running a cell optimization or NpT molecular dynamics. "// & "Check the &CELL section for further details.") CALL create_cell_section_low(subsection, periodic) CALL section_add_subsection(section, subsection) diff --git a/src/mscfg_methods.F b/src/mscfg_methods.F index dd87e58696..88721beab6 100644 --- a/src/mscfg_methods.F +++ b/src/mscfg_methods.F @@ -18,7 +18,6 @@ MODULE mscfg_methods USE almo_scf_types, ONLY: almo_scf_env_type USE atomic_kind_types, ONLY: get_atomic_kind - USE cell_types, ONLY: cell_type USE cp_dbcsr_api, ONLY: dbcsr_copy,& dbcsr_create,& dbcsr_type_no_symmetry @@ -151,7 +150,6 @@ CONTAINS INTEGER, DIMENSION(:), POINTER :: atom_index, atom_list LOGICAL :: global_almo_scf_keyword, smear_almo_scf TYPE(almo_scf_env_type), POINTER :: almo_scf_env - TYPE(cell_type), POINTER :: cell TYPE(cp_subsys_type), POINTER :: subsys, subsys_loc TYPE(mo_set_type), DIMENSION(:), POINTER :: mos, mos_of_frag TYPE(molecular_scf_guess_env_type), POINTER :: mscfg_env @@ -165,7 +163,7 @@ CONTAINS CALL timeset(routineN, handle) - NULLIFY (subsys_loc, subsys, particles, para_env, cell, atom_index, atom_type, & + NULLIFY (subsys_loc, subsys, particles, para_env, atom_index, atom_type, & force_env_section, qs_env_loc, mscfg_env, qs_env, qs_energy) CALL force_env_get(force_env, force_env_section=force_env_section, & qs_env=qs_env) @@ -185,8 +183,7 @@ CONTAINS ! ! Get access to critical data before the loop ! - CALL force_env_get(force_env=force_env, subsys=subsys, para_env=para_env, & - cell=cell) + CALL force_env_get(force_env=force_env, subsys=subsys, para_env=para_env) CALL cp_subsys_get(subsys, particles=particles) CALL get_qs_env(qs_env, mscfg_env=mscfg_env, almo_scf_env=almo_scf_env) CPASSERT(ASSOCIATED(mscfg_env)) @@ -280,8 +277,8 @@ CONTAINS ! ! Create the environment of a subsystem ! - CALL create_small_subsys(subsys_loc, big_subsys=subsys, & - small_para_env=para_env, small_cell=cell, sub_atom_index=atom_index, & + CALL create_small_subsys(subsys_loc, big_subsys=subsys, small_para_env=para_env, & + small_cell=subsys%cell, sub_atom_index=atom_index, & sub_atom_kind_name=atom_type, para_env=para_env, & force_env_section=force_env_section, subsys_section=subsys_section) ALLOCATE (qs_env_loc) diff --git a/src/qmmm_create.F b/src/qmmm_create.F index c8857ca1ca..ae43a5ef42 100644 --- a/src/qmmm_create.F +++ b/src/qmmm_create.F @@ -16,7 +16,8 @@ MODULE qmmm_create Laino2005,& Laino2006,& cite_reference - USE cell_methods, ONLY: write_cell + USE cell_methods, ONLY: read_cell,& + write_cell USE cell_types, ONLY: cell_clone,& cell_release,& cell_type,& @@ -141,13 +142,13 @@ CONTAINS TYPE(qmmm_env_qm_type), POINTER :: qmmm_env_qm TYPE(qmmm_links_type), POINTER :: qmmm_links TYPE(qs_environment_type), POINTER :: qs_env - TYPE(section_vals_type), POINTER :: multipole_section, print_gen, & - print_section, qmmm_periodic + TYPE(section_vals_type), POINTER :: cell_section, multipole_section, & + print_gen, print_section, qmmm_periodic CALL timeset(routineN, handle) NULLIFY (qm_atom_index, mm_atom_index, qm_atom_type) - NULLIFY (qmmm_env_qm, subsys_mm, subsys_qm, mm_cell, qm_cell_small) + NULLIFY (qmmm_env_qm, subsys_mm, subsys_qm, mm_cell, qm_cell_small, cell_section) NULLIFY (mm_atom_chrg, mm_el_pot_radius, qmmm_env_mm, fist_env) NULLIFY (mm_link_atoms, mm_link_scale_factor, qmmm_links, added_charges, added_shells) NULLIFY (fist_scale_charge_link, print_section, fist_nonbond_env) @@ -242,8 +243,15 @@ CONTAINS CALL fist_env_set(fist_env, qmmm=.TRUE., qmmm_env=qmmm_env_mm) CALL fist_init(fist_env, root_section, para_env, force_env_section, & subsys_section, use_motion_section, prev_subsys=prev_subsys) + + ! Get cell: mm_cell from FORCE_EVAL/SUBSYS/CELL like in MM-only Fist; + ! qm_cell from FORCE_EVAL/QMMM/CELL CALL fist_env_get(fist_env, subsys=subsys_mm, cell=mm_cell) mm_cell%tag = "CELL_MM" + cell_section => section_vals_get_subs_vals(qmmm_section, "CELL") + CALL read_cell(qm_cell_small, qm_cell_small, cell_section=cell_section, & + check_for_ref=.FALSE., para_env=para_env) + qm_cell_small%tag = "CELL_QM" ! Set up QM/MM Options CALL setup_qmmm_vars_qm(qmmm_section, & @@ -360,8 +368,9 @@ CONTAINS ALLOCATE (qs_env) CALL qs_env_create(qs_env, globenv) CALL qs_init(qs_env, para_env, root_section, globenv=globenv, cp_subsys=subsys_qm, & - cell=qm_cell_small, qmmm=.TRUE., qmmm_env_qm=qmmm_env_qm, & - force_env_section=force_env_section, subsys_section=subsys_section, & + qmmm=.TRUE., qmmm_env_qm=qmmm_env_qm, & + force_env_section=force_env_section, & + subsys_section=subsys_section, & use_motion_section=use_motion_section) CALL cp_subsys_release(subsys_qm) diff --git a/src/qmmm_init.F b/src/qmmm_init.F index 8e8ba81e00..c9accbb2e2 100644 --- a/src/qmmm_init.F +++ b/src/qmmm_init.F @@ -16,7 +16,6 @@ MODULE qmmm_init USE atomic_kind_types, ONLY: atomic_kind_type,& get_atomic_kind,& set_atomic_kind - USE cell_methods, ONLY: read_cell USE cell_types, ONLY: cell_type,& use_perd_xyz USE cp_control_types, ONLY: dft_control_type @@ -649,22 +648,15 @@ CONTAINS TYPE(atomic_kind_list_type), POINTER :: atomic_kinds TYPE(atomic_kind_type), POINTER :: atomic_kind TYPE(fist_potential_type), POINTER :: fist_potential - TYPE(section_vals_type), POINTER :: cell_section, eri_mme_section, & - image_charge_section, mm_kinds + TYPE(section_vals_type), POINTER :: eri_mme_section, image_charge_section, & + mm_kinds - NULLIFY (mm_link_atoms, cell_section, tmp_sph_cut) + NULLIFY (mm_link_atoms, tmp_sph_cut) NULLIFY (image_charge_section) qmmm_link = .FALSE. CALL section_vals_get(qmmm_section, explicit=explicit) IF (explicit) THEN - ! - ! QM_CELL - ! - cell_section => section_vals_get_subs_vals(qmmm_section, "CELL") - CALL read_cell(qm_cell_small, qm_cell_small, cell_section=cell_section, & - check_for_ref=.FALSE., para_env=para_env) - qm_cell_small%tag = "CELL_QM" CALL section_vals_val_get(qmmm_section, "E_COUPL", i_val=qmmm_coupl_type) CALL section_vals_val_get(qmmm_section, "EPS_MM_RSPACE", r_val=eps_mm_rspace) CALL section_vals_val_get(qmmm_section, "SPHERICAL_CUTOFF", r_vals=tmp_sph_cut) diff --git a/src/qs_environment.F b/src/qs_environment.F index 9aa8d23736..c01b436806 100644 --- a/src/qs_environment.F +++ b/src/qs_environment.F @@ -240,8 +240,6 @@ CONTAINS !> \param globenv ... !> \param cp_subsys ... !> \param kpoint_env ... -!> \param cell ... -!> \param cell_ref ... !> \param qmmm ... !> \param qmmm_env_qm ... !> \param force_env_section ... @@ -250,7 +248,7 @@ CONTAINS !> \param silent ... !> \author Creation (22.05.2000,MK) ! ************************************************************************************************** - SUBROUTINE qs_init(qs_env, para_env, root_section, globenv, cp_subsys, kpoint_env, cell, cell_ref, & + SUBROUTINE qs_init(qs_env, para_env, root_section, globenv, cp_subsys, kpoint_env, & qmmm, qmmm_env_qm, force_env_section, subsys_section, & use_motion_section, silent) @@ -260,7 +258,6 @@ CONTAINS TYPE(global_environment_type), OPTIONAL, POINTER :: globenv TYPE(cp_subsys_type), OPTIONAL, POINTER :: cp_subsys TYPE(kpoint_type), OPTIONAL, POINTER :: kpoint_env - TYPE(cell_type), OPTIONAL, POINTER :: cell, cell_ref LOGICAL, INTENT(IN), OPTIONAL :: qmmm TYPE(qmmm_env_qm_type), OPTIONAL, POINTER :: qmmm_env_qm TYPE(section_vals_type), POINTER :: force_env_section, subsys_section @@ -336,7 +333,7 @@ CONTAINS subsys_section=subsys_section, & use_motion_section=use_motion_section, & root_section=root_section, & - cp_subsys=cp_subsys, cell=cell, cell_ref=cell_ref, & + cp_subsys=cp_subsys, & elkind=is_semi, silent=silent) ALLOCATE (ks_env) diff --git a/src/qs_subsys_methods.F b/src/qs_subsys_methods.F index e16139b937..fc352effb5 100644 --- a/src/qs_subsys_methods.F +++ b/src/qs_subsys_methods.F @@ -15,16 +15,9 @@ MODULE qs_subsys_methods get_atomic_kind USE basis_set_types, ONLY: get_gto_basis_set,& gto_basis_set_type - USE cell_methods, ONLY: cell_create,& - read_cell,& - write_cell - USE cell_types, ONLY: cell_clone,& - cell_release,& - cell_type USE cp_subsys_methods, ONLY: cp_subsys_create USE cp_subsys_types, ONLY: cp_subsys_get,& cp_subsys_release,& - cp_subsys_set,& cp_subsys_type USE external_potential_types, ONLY: all_potential_type,& get_potential,& @@ -65,30 +58,26 @@ CONTAINS !> \param subsys_section ... !> \param use_motion_section ... !> \param cp_subsys ... -!> \param cell ... -!> \param cell_ref ... !> \param elkind ... !> \param silent ... ! ************************************************************************************************** SUBROUTINE qs_subsys_create(subsys, para_env, root_section, force_env_section, subsys_section, & - use_motion_section, cp_subsys, cell, cell_ref, elkind, silent) + use_motion_section, cp_subsys, elkind, silent) TYPE(qs_subsys_type), INTENT(OUT) :: subsys TYPE(mp_para_env_type), POINTER :: para_env TYPE(section_vals_type), OPTIONAL, POINTER :: root_section TYPE(section_vals_type), POINTER :: force_env_section, subsys_section LOGICAL, INTENT(IN) :: use_motion_section TYPE(cp_subsys_type), OPTIONAL, POINTER :: cp_subsys - TYPE(cell_type), OPTIONAL, POINTER :: cell, cell_ref LOGICAL, INTENT(IN), OPTIONAL :: elkind, silent - LOGICAL :: be_silent, use_ref_cell + LOGICAL :: be_silent TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set - TYPE(cell_type), POINTER :: my_cell, my_cell_ref TYPE(cp_subsys_type), POINTER :: my_cp_subsys TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set - TYPE(section_vals_type), POINTER :: cell_section, kind_section + TYPE(section_vals_type), POINTER :: kind_section - NULLIFY (atomic_kind_set, qs_kind_set, cell_section, kind_section, my_cell, my_cell_ref, my_cp_subsys) + NULLIFY (atomic_kind_set, qs_kind_set, kind_section, my_cp_subsys) be_silent = .FALSE. IF (PRESENT(silent)) be_silent = silent @@ -105,37 +94,6 @@ CONTAINS CPABORT("qs_subsys_create: cp_subsys or root_section needed") END IF - ! create cp_subsys%cell - !TODO: moved to cp_subsys_create(), needs further disentanglement of cell_ref. - ![NB] subroutines qs_init() and qs_init_subsys() may pass cell and cell_ref to - ! qs_subsys_create() here, and two places where qs_init() is called are tricky: - ! 1. For QM/MM calculation, qmmm_env_create() passes "qm_cell_small" from - ! setup_qmmm_vars_qm() to qs_init() and subsequently to my_cell here; - ! 2. For BSSE calculation, eval_bsse_energy_low() calls create_small_subsys() - ! to create cp_subsys_type "subsys_loc" and then calls qs_init() with - ! subsys_loc, but create_small_subsys() does not set subsys_loc%cell. - ! If not for these, it would have been easy to set my_cell => my_cp_subsys%cell - ! and my_cell_ref => my_cp_subsys%cell_ref so as to avoid re-read the cell_section - ! when cell and cell_ref are not passed as input here. - use_ref_cell = .FALSE. - IF (PRESENT(cell)) THEN - my_cell => cell - IF (PRESENT(cell_ref)) THEN - my_cell_ref => cell_ref - use_ref_cell = .TRUE. - ELSE - CALL cell_create(my_cell_ref) - CALL cell_clone(my_cell, my_cell_ref, tag="CELL_REF") - END IF - ELSE - cell_section => section_vals_get_subs_vals(subsys_section, "CELL") - CALL read_cell(my_cell, my_cell_ref, use_ref_cell=use_ref_cell, & - cell_section=cell_section, para_env=para_env) - END IF - CALL cp_subsys_set(my_cp_subsys, cell=my_cell) - CALL write_cell(my_cell, subsys_section) - CALL write_cell(my_cell_ref, subsys_section) - ! setup qs_kinds CALL cp_subsys_get(my_cp_subsys, atomic_kind_set=atomic_kind_set) kind_section => section_vals_get_subs_vals(subsys_section, "KIND") @@ -147,12 +105,10 @@ CONTAINS CALL qs_subsys_set(subsys, & cp_subsys=my_cp_subsys, & - cell_ref=my_cell_ref, & - use_ref_cell=use_ref_cell, & + cell_ref=my_cp_subsys%cell_ref, & + use_ref_cell=my_cp_subsys%use_ref_cell, & qs_kind_set=qs_kind_set) - IF (.NOT. PRESENT(cell)) CALL cell_release(my_cell) - IF (.NOT. PRESENT(cell_ref)) CALL cell_release(my_cell_ref) IF (.NOT. PRESENT(cp_subsys)) CALL cp_subsys_release(my_cp_subsys) END SUBROUTINE qs_subsys_create diff --git a/src/topology.F b/src/topology.F index 91bb37ff77..a6002ed3ba 100644 --- a/src/topology.F +++ b/src/topology.F @@ -152,17 +152,17 @@ CONTAINS INTEGER :: handle, iw, iw2 LOGICAL :: binary_coord_read, el_as_kind, explicit, & - my_qmmm + my_qmmm, use_ref_cell TYPE(cell_type), POINTER :: my_cell, my_cell_ref TYPE(cp_logger_type), POINTER :: logger - TYPE(section_vals_type), POINTER :: cell_section, constraint_section, & - topology_section + TYPE(section_vals_type), POINTER :: cell_ref_section, cell_section, & + constraint_section, topology_section TYPE(topology_parameters_type) :: topology NULLIFY (logger) logger => cp_get_default_logger() CALL timeset(routineN, handle) - NULLIFY (cell_section, constraint_section, topology_section) + NULLIFY (cell_section, cell_ref_section, constraint_section, topology_section) cell_section => section_vals_get_subs_vals(subsys_section, "CELL") IF (use_motion_section) THEN @@ -192,13 +192,17 @@ CONTAINS NULLIFY (my_cell, my_cell_ref) CALL read_cell(my_cell, my_cell_ref, cell_section=cell_section, & topology_section=topology_section, para_env=para_env) + use_ref_cell = .FALSE. + cell_ref_section => section_vals_get_subs_vals(cell_section, "CELL_REF") + CALL section_vals_get(cell_ref_section, explicit=use_ref_cell) IF (.NOT. ASSOCIATED(topology%cell)) CALL cell_create(topology%cell) IF (.NOT. ASSOCIATED(topology%cell_ref)) CALL cell_create(topology%cell_ref) CALL cell_clone(my_cell, topology%cell) CALL cell_clone(my_cell_ref, topology%cell_ref) CALL setup_cell_muc(topology%cell_muc, topology%cell, subsys_section) IF (PRESENT(subsys)) THEN ! subsys is passed from cp_subsys_create - CALL cp_subsys_set(subsys, cell=my_cell, cell_ref=my_cell_ref) + CALL cp_subsys_set(subsys, cell=my_cell, cell_ref=my_cell_ref, & + use_ref_cell=use_ref_cell) END IF CALL write_cell(my_cell, subsys_section, tag="CELL") CALL write_cell(my_cell_ref, subsys_section, tag="CELL_REF")