From 70a54851ac1fbccc9fbe91c3958dbac8ca71dd32 Mon Sep 17 00:00:00 2001 From: Teodoro Laino Date: Mon, 26 Mar 2012 09:02:58 +0000 Subject: [PATCH] Preliminary framework for adaptive QM/MM (Noam and Teo) + 1 regtest svn-origin-rev: 12203 --- src/OBJECTDEFS | 1 + src/cell_types.F | 23 +- src/force_env_methods.F | 103 ++- src/force_env_types.F | 31 +- src/input_cp2k_qmmm.F | 176 ++++- src/input_cp2k_restarts.F | 3 +- src/integrator.F | 30 + src/qmmm_force_mixing.F | 735 ++++++++++++++++++ src/qmmm_main.F | 210 ++++- src/qmmm_types.F | 9 + src/qmmm_util.F | 139 +++- .../H2O-qmmm-none-force-mixing-1.inp | 122 +++ tests/QMMM/QS/regtest-1/TEST_FILES | 1 + tests/QMMM/QS/sample_pdb/5H2O.pdb | 20 + tests/QMMM/QS/sample_pot/5H2O.pot | 76 ++ tests/QMMM/QS/sample_psf/5H2O.psf | 45 ++ tests/QMMM/QS/sample_xyz/5H2O_adaptive.xyz | 136 ++++ 17 files changed, 1813 insertions(+), 47 deletions(-) create mode 100644 src/qmmm_force_mixing.F create mode 100644 tests/QMMM/QS/regtest-1/H2O-qmmm-none-force-mixing-1.inp create mode 100644 tests/QMMM/QS/sample_pdb/5H2O.pdb create mode 100644 tests/QMMM/QS/sample_pot/5H2O.pot create mode 100644 tests/QMMM/QS/sample_psf/5H2O.psf create mode 100644 tests/QMMM/QS/sample_xyz/5H2O_adaptive.xyz diff --git a/src/OBJECTDEFS b/src/OBJECTDEFS index 91db538bc4..a86c9ef1b4 100644 --- a/src/OBJECTDEFS +++ b/src/OBJECTDEFS @@ -431,6 +431,7 @@ OBJECTS_GENERIC =\ pw_types.o\ qmmm_elpot.o\ qmmm_ff_fist.o\ + qmmm_force_mixing.o\ qmmm_gaussian_data.o\ qmmm_gaussian_init.o\ qmmm_gaussian_input.o\ diff --git a/src/cell_types.F b/src/cell_types.F index fe7ff2ac12..052d4f48f3 100644 --- a/src/cell_types.F +++ b/src/cell_types.F @@ -68,7 +68,7 @@ MODULE cell_types ! *** Public subroutines *** PUBLIC :: get_cell, get_cell_param, init_cell, read_cell,& write_cell, cell_create, cell_retain, cell_release,& - cell_clone, compare_cells, parse_cell_line, set_cell_param + cell_clone, cell_copy, compare_cells, parse_cell_line, set_cell_param ! *** Public functions *** PUBLIC :: plane_distance, pbc, real_to_scaled, scaled_to_real @@ -131,6 +131,27 @@ CONTAINS cell_out%id_nr = last_cell_id END SUBROUTINE cell_clone +! ***************************************************************************** + SUBROUTINE cell_copy (cell_in, cell_out, error) + TYPE(cell_type), POINTER :: cell_in, cell_out + TYPE(cp_error_type), INTENT(INOUT) :: error + + CHARACTER(len=*), PARAMETER :: routineN = 'cell_copy', & + routineP = moduleN//':'//routineN + + LOGICAL :: failure + + failure=.FALSE. + CPPrecondition(ASSOCIATED(cell_in),cp_failure_level,routineP,error,failure) + CPPrecondition(ASSOCIATED(cell_out),cp_failure_level,routineP,error,failure) + + cell_out%deth = cell_in%deth + cell_out%perd = cell_in%perd + cell_out%hmat = cell_in%hmat + cell_out%h_inv = cell_in%h_inv + cell_out%orthorhombic = cell_in%orthorhombic + END SUBROUTINE cell_copy + ! ***************************************************************************** !> \brief Read cell info from a line (parsed from a file) !> \author Teodoro Laino [tlaino] - University of Zurich diff --git a/src/force_env_methods.F b/src/force_env_methods.F index ce0b342526..21f977f1c4 100644 --- a/src/force_env_methods.F +++ b/src/force_env_methods.F @@ -112,11 +112,10 @@ MODULE force_env_methods qmmm_added_chrg_forces,& qmmm_link_Imomm_coord,& qmmm_link_Imomm_forces - USE qmmm_types, ONLY: fist_subsys,& - qmmm_env_qm_retain,& - qmmm_env_qm_type,& - qmmm_links_type,& - qs_subsys + USE qmmm_types, ONLY: & + fist_subsys, force_mixing_core_subsys, force_mixing_extended_subsys, & + force_mixing_label_QM_extended, primary_subsys, qmmm_env_qm_retain, & + qmmm_env_qm_type, qmmm_links_type, qs_subsys USE qmmm_util, ONLY: apply_qmmm_translate,& apply_qmmm_walls USE qs_energy, ONLY: qs_energies @@ -677,7 +676,7 @@ CONTAINS force_env%in_use=use_qmmm force_env%qmmm_env => qmmm_env CALL qmmm_env_qm_retain(qmmm_env,error=error) - force_env%virial => sub_force_env(1)%force_env%virial + force_env%virial => sub_force_env(primary_subsys)%force_env%virial CALL virial_retain(force_env%virial,error=error) ! Virial controlled through the external request CALL virial_set ( virial=force_env%virial,& @@ -734,7 +733,7 @@ CONTAINS CALL eip_env_get(force_env%eip_env,subsys=force_env%subsys,error=error) CALL cp_subsys_retain(force_env%subsys, error=error) CASE(use_qmmm) - subsys => force_env%sub_force_env(1)%force_env%subsys + subsys => force_env%sub_force_env(primary_subsys)%force_env%subsys force_env%subsys => subsys CALL cp_subsys_retain(subsys,error=error) CASE(use_mixed_force) @@ -794,7 +793,6 @@ CONTAINS !> \author Fawzi Mohamed ! ***************************************************************************** RECURSIVE SUBROUTINE qmmm_energy_and_forces(force_env,calc_force,error) - TYPE(force_env_type), POINTER :: force_env LOGICAL, INTENT(IN), OPTIONAL :: calc_force TYPE(cp_error_type), INTENT(inout) :: error @@ -802,6 +800,90 @@ CONTAINS CHARACTER(len=*), PARAMETER :: routineN = 'qmmm_energy_and_forces', & routineP = moduleN//':'//routineN + INTEGER :: ip, isubf + INTEGER, POINTER :: cur_indices(:), cur_labels(:) + LOGICAL :: failure + TYPE(cp_subsys_type), POINTER :: subsys_primary, & + subsys_qmmm_core, & + subsys_qmmm_extended + TYPE(particle_type), DIMENSION(:), & + POINTER :: particles_primary, & + particles_qmmm_core, & + particles_qmmm_extended + TYPE(section_vals_type), POINTER :: force_env_section + +! Possibly translate the system + + CALL apply_qmmm_translate(force_env, error) + + DO isubf=1, SIZE(force_env%sub_force_env) + CALL qmmm_energy_and_forces_low(force_env%sub_force_env(isubf)%force_env,& + calc_force,error) + END DO + + IF (SIZE(force_env%sub_force_env) == 1) THEN + ! see force_env_get for potential energy to copy + ! copy forces, energy, virial up + CONTINUE + ELSE IF (SIZE(force_env%sub_force_env) == 2) THEN + ! force mixing happens here + ! energy will just inherit from extended sub_force_env in force_env_get() + + ! get forces from subsys of each sub force env + CALL force_env_get(force_env%sub_force_env(force_mixing_core_subsys)%force_env,& + subsys=subsys_qmmm_core,force_env_section=force_env_section,error=error) + CALL force_env_get(force_env%sub_force_env(force_mixing_extended_subsys)%force_env,& + subsys=subsys_qmmm_extended,error=error) + + CALL section_vals_val_get(force_env_section,"QMMM%FORCE_MIXING%RESTART_INFO%INDICES",i_vals=cur_indices,error=error) + CALL section_vals_val_get(force_env_section,"QMMM%FORCE_MIXING%RESTART_INFO%LABELS",i_vals=cur_labels,error=error) + + particles_qmmm_extended => subsys_qmmm_extended%particles%els + particles_qmmm_core => subsys_qmmm_core%particles%els + DO ip=1,SIZE(cur_indices) + IF (cur_labels(ip) >= force_mixing_label_QM_extended) THEN ! this is a QM atom + ! copy (QM) force from extended calculation + particles_qmmm_core(cur_indices(ip))%f=particles_qmmm_extended(cur_indices(ip))%f + ENDIF + END DO + + IF (force_mixing_core_subsys /= primary_subsys) THEN + CALL force_env_get(force_env%sub_force_env(primary_subsys)%force_env,& + subsys=subsys_primary,force_env_section=force_env_section,error=error) + particles_primary => subsys_primary%particles%els + DO ip=1,SIZE(particles_qmmm_core) + particles_primary(ip)%f=particles_qmmm_core(ip)%f + END DO + ENDIF + + ELSE + CALL cp_assert(.FALSE.,cp_failure_level,cp_assertion_failed,& + routineP,"size(sub_force_env) neither 1 (conventional QM/MM) nor "//& + "2 (force-mixing)! "//& +CPSourceFileRef,& + error,failure) + ENDIF + + END SUBROUTINE qmmm_energy_and_forces + +! ***************************************************************************** +!> \brief calculates the qm/mm energy and forces +!> \param calc_force if also the forces should be calculated +!> \param error variable to control error logging, stopping,... +!> see module cp_error_handling +!> \par History +!> 05.2004 created [fawzi] +!> \author Fawzi Mohamed +! ***************************************************************************** + RECURSIVE SUBROUTINE qmmm_energy_and_forces_low(force_env,calc_force,error) + + TYPE(force_env_type), POINTER :: force_env + LOGICAL, INTENT(IN), OPTIONAL :: calc_force + TYPE(cp_error_type), INTENT(inout) :: error + + CHARACTER(len=*), PARAMETER :: routineN = 'qmmm_energy_and_forces_low', & + routineP = moduleN//':'//routineN + INTEGER :: ip, output_unit INTEGER, DIMENSION(:), POINTER :: qm_atom_index LOGICAL :: calculate_forces, failure, & @@ -849,9 +931,6 @@ CONTAINS END IF CPPrecondition(ASSOCIATED(qm_atom_index),cp_failure_level,routineP,error,failure) - ! Possibly translate the system - CALL apply_qmmm_translate(force_env, error) - particles_mm => subsys_mm%particles%els particles_qm => subsys_qm%particles%els @@ -925,7 +1004,7 @@ CONTAINS END IF CALL cp_print_key_finished_output(output_unit,logger,force_env_section,& "QMMM%PRINT%DERIVATIVES",error=error) - END SUBROUTINE qmmm_energy_and_forces + END SUBROUTINE qmmm_energy_and_forces_low ! ***************************************************************************** !> \brief ****f* force_env_methods/mixed_energy_forces [1.0] diff --git a/src/force_env_types.F b/src/force_env_types.F index 34a5f7690f..96d1693b03 100644 --- a/src/force_env_types.F +++ b/src/force_env_types.F @@ -60,8 +60,12 @@ MODULE force_env_types USE mixed_environment_types, ONLY: get_mixed_env,& mixed_env_release,& mixed_environment_type - USE qmmm_types, ONLY: qmmm_env_qm_release,& - qmmm_env_qm_type + USE qmmm_types, ONLY: fist_subsys,& + force_mixing_extended_subsys,& + primary_subsys,& + qmmm_env_qm_release,& + qmmm_env_qm_type,& + qs_subsys USE qs_energy_types, ONLY: qs_energy_type USE qs_environment_methods, ONLY: qs_env_rebuild_pw_env USE qs_environment_types, ONLY: get_qs_env,& @@ -355,6 +359,7 @@ CONTAINS CHARACTER(len=*), PARAMETER :: routineN = 'force_env_get', & routineP = moduleN//':'//routineN + INTEGER :: cur_subsys LOGICAL :: failure REAL(KIND=dp) :: eip_kinetic_energy, & eip_potential_energy, & @@ -426,26 +431,38 @@ CONTAINS CPPrecondition(.NOT.PRESENT(kinetic_energy),cp_failure_level,routineP,error,failure) CASE (use_qmmm) IF (PRESENT(cell)) THEN - CALL force_env_get(force_env%sub_force_env(1)%force_env,& + CALL force_env_get(force_env%sub_force_env(primary_subsys)%force_env,& cell=cell,& error=error) ENDIF IF (PRESENT(cell_ref)) THEN - CALL force_env_get(force_env%sub_force_env(1)%force_env,& + CALL force_env_get(force_env%sub_force_env(primary_subsys)%force_env,& cell_ref=cell_ref,& error=error) ENDIF IF (PRESENT(use_ref_cell)) use_ref_cell = .FALSE. IF (PRESENT(kinetic_energy)) THEN - CALL force_env_get(force_env%sub_force_env(1)%force_env,& + CALL force_env_get(force_env%sub_force_env(primary_subsys)%force_env,& kinetic_energy=kinetic_energy,& error=error) END IF IF (PRESENT(potential_energy)) THEN - CALL force_env_get(force_env%sub_force_env(1)%force_env,& + ! get the underlying energies from primary subsys. This is the only subsys + ! for conventional QM/MM, and force-mixing knows to put relevant energy there. + IF (SIZE(force_env%sub_force_env) == 1) THEN + cur_subsys = primary_subsys + ELSE IF (SIZE(force_env%sub_force_env) == 2) THEN + cur_subsys = force_mixing_extended_subsys + ELSE + CALL cp_assert(.FALSE.,cp_failure_level,cp_assertion_failed,& + routineP,"bad number of sub_force_envs getting qmmm potential energy "//& + CPSourceFileRef,& + error,failure) + ENDIF + CALL force_env_get(force_env%sub_force_env(cur_subsys)%force_env%sub_force_env(fist_subsys)%force_env,& potential_energy=penergy_mm,& error=error) - CALL force_env_get(force_env%sub_force_env(2)%force_env,& + CALL force_env_get(force_env%sub_force_env(cur_subsys)%force_env%sub_force_env(qs_subsys)%force_env,& potential_energy=penergy_qm,& error=error) potential_energy = penergy_qm+penergy_mm diff --git a/src/input_cp2k_qmmm.F b/src/input_cp2k_qmmm.F index 65c627812d..9bec8a2965 100644 --- a/src/input_cp2k_qmmm.F +++ b/src/input_cp2k_qmmm.F @@ -36,7 +36,8 @@ MODULE input_cp2k_qmmm section_type USE input_val_types, ONLY: char_t,& integer_t,& - real_t + real_t, & + logical_t USE kinds, ONLY: dp USE string_utilities, ONLY: s2a #include "cp_common_uses.h" @@ -189,6 +190,11 @@ CONTAINS CALL section_add_keyword(section,keyword,error=error) CALL keyword_release(keyword,error=error) + ! NB: remember to create these + CALL create_qmmm_force_mixing_section(subsection, error) + CALL section_add_subsection(section,subsection,error=error) + CALL section_release(subsection,error=error) + CALL create_qmmm_qm_kinds(subsection,error) CALL section_add_subsection(section,subsection,error=error) CALL section_release(subsection,error=error) @@ -278,6 +284,174 @@ CONTAINS END SUBROUTINE create_qmmm_mm_kinds +! ***************************************************************************** +!> \brief Input section to create FORCE_MIXING sections +!> \param error variable to control error logging, stopping,... +!> see module cp_error_handling +!> \author noam +! ***************************************************************************** + SUBROUTINE create_qmmm_force_mixing_section(section, error) + TYPE(section_type), POINTER :: section + TYPE(cp_error_type), INTENT(inout) :: error + + CHARACTER(len=*), PARAMETER :: & + routineN = 'create_qmmm_force_mixing_section', & + routineP = moduleN//':'//routineN + + LOGICAL :: failure + TYPE(keyword_type), POINTER :: keyword + TYPE(section_type), POINTER :: link_subsection, print_key, & + subsection + + failure=.FALSE. + NULLIFY(keyword) + CPPrecondition(.NOT.ASSOCIATED(section),cp_failure_level,routineP,error,failure) + CALL section_create(section,name="FORCE_MIXING",& + description="This section defines parameters for force-mixing based QM/MM,"//& + " which actually does two conventional QM/MM calculations, on a small "//& + " and a large QM regions, and combines the MM forces from one and QM "//& + " forces from the other to create a complete set of forces. Energy is "//& + " not conserved (although the energy from the large QM region is saved) "//& + " so a proper thermostat (i.e. massive, and able to handle dissipation) must be used.",& + n_keywords=5, n_subsections=3, repeats=.FALSE., required=.FALSE.,& + error=error) + + CALL keyword_create(keyword, name="_SECTION_PARAMETERS_", & + description="Enables force-mixing", & + default_l_val=.FALSE., lone_keyword_l_val=.TRUE., error=error) + CALL section_add_keyword(section,keyword,error=error) + CALL keyword_release(keyword,error=error) + + CALL keyword_create(keyword, name="R_CORE",& + description="Specify the inner and outer radii of core QM region. "//& + "All atoms within this distance (hysteretically) of any atoms listed in "//& + "ATOM_INDEX will be core QM atoms in the force-mixing calculation.", & + usage="R_CORE ",n_var=2,type_of_var=real_t,& + default_r_vals=(/ cp_unit_to_cp2k(0.0_dp,"angstrom",error=error),& + cp_unit_to_cp2k(0.0_dp,"angstrom",error=error) /),& + unit_str="angstrom",error=error) + CALL section_add_keyword(section,keyword,error=error) + CALL keyword_release(keyword,error=error) + + CALL keyword_create(keyword, name="R_QM",& + description="Specify the inner and outer radii of extended QM region. "//& + "All atoms within this distance (hysteretically) of any atoms listed in "//& + "ATOM_INDEX will be QM atoms in the force-mixing calculation.", & + usage="R_QM ",n_var=2,type_of_var=real_t,& + default_r_vals=(/ cp_unit_to_cp2k(0.5_dp,"angstrom",error=error),& + cp_unit_to_cp2k(1.0_dp,"angstrom",error=error) /),& + unit_str="angstrom",error=error) + CALL section_add_keyword(section,keyword,error=error) + CALL keyword_release(keyword,error=error) + + CALL keyword_create(keyword, name="QM_EXTENDED_SEED_IS_ONLY_CORE_LIST",& + description="Makes the extended QM zone be defined hysterestically "//& + "by distance from QM core list (i.e. atoms specified explicitly by "//& + "user) instead of from full QM core region (specified by user + hysteretic "//& + "selection + unbreakable bonds)",& + usage="QM_EXTENDED_SEED_IS_ONLY_CORE_LIST ",n_var=1,type_of_var=logical_t,& + default_l_val=.FALSE., repeats = .FALSE., required = .FALSE., error=error) + CALL section_add_keyword(section,keyword,error=error) + CALL keyword_release(keyword,error=error) + + CALL keyword_create(keyword, name="R_BUF",& + description="Specify the inner and outer radii of buffer region. "//& + "All atoms within this distance (hysteretically) of any QM atoms "//& + "will be buffer atoms in the force-mixing calculation.", & + usage="R_BUF ",n_var=2,type_of_var=real_t,& + default_r_vals=(/ cp_unit_to_cp2k(0.5_dp,"angstrom",error=error),& + cp_unit_to_cp2k(1.0_dp,"angstrom",error=error) /),& + unit_str="angstrom",error=error) + CALL section_add_keyword(section,keyword,error=error) + CALL keyword_release(keyword,error=error) + + CALL keyword_create(keyword, name="QM_KIND_ELEMENT_MAPPING",& + description="Mapping from elements to QM_KINDs for adaptivity",& + usage="QM_KIND_ELEMENT_MAPPING {El} {QM_KIND}",& + n_var=2,type_of_var=char_t,repeats=.TRUE.,required=.TRUE.,error=error) + CALL section_add_keyword(section,keyword,error=error) + CALL keyword_release(keyword,error=error) + + CALL keyword_create(keyword, name="MAX_N_QM",& + description="Maximum number of QM atoms, for detection of runaway adaptive selection.",& + usage="MAX_N_QM int",default_i_val=300,& + n_var=1,type_of_var=integer_t,repeats=.FALSE.,required=.FALSE.,error=error) + CALL section_add_keyword(section,keyword,error=error) + CALL keyword_release(keyword,error=error) + + + ![NB] also need a list? + ![NB] maybe not list+links , but some sort of link template + ![NB] also, breakable bonds? + ! BUFFER_LINKS subsection + NULLIFY(subsection) + CALL section_create(subsection,name="BUFFER_LINKS",& + description="Information about possible links for the buffer QM/MM calculation."//& + "Not used yet - need to implement buffer selection by atom and walking of connectivity data.",& + n_keywords=0, n_subsections=1, repeats=.TRUE., required=.TRUE.,& + error=error) + + NULLIFY(link_subsection) + CALL create_qmmm_link_section(link_subsection,error) + CALL section_add_subsection(subsection,link_subsection,error=error) + CALL section_release(link_subsection,error=error) + + CALL section_add_subsection(section,subsection,error=error) + CALL section_release(subsection,error=error) + + ! RESTART_INFO subsection + NULLIFY(subsection) + CALL section_create(subsection,name="RESTART_INFO",& + description="This section provides information about old force-mixing indices and labels, "//& + "for restarts.", & + n_keywords=2, n_subsections=0, repeats=.FALSE., required=.TRUE.,& + error=error) + + CALL keyword_create(keyword, name="INDICES",& + description="Indices of atoms in previous step QM regions.", & + usage="INDICES 1 2 ...", & + n_var=-1,type_of_var=integer_t,repeats=.TRUE.,required=.TRUE.,error=error) + CALL section_add_keyword(subsection,keyword,error=error) + CALL keyword_release(keyword,error=error) + + CALL keyword_create(keyword, name="LABELS",& + description="Labels of atoms in previous step QM regions.", & + usage="LABELS 1 1 ...", & + n_var=-1,type_of_var=integer_t,repeats=.TRUE.,required=.TRUE.,error=error) + CALL section_add_keyword(subsection,keyword,error=error) + CALL keyword_release(keyword,error=error) + + CALL section_add_subsection(section,subsection,error=error) + CALL section_release(subsection,error=error) + + ! PRINT subsection, with keys for neighbor list + CALL section_create(subsection,name="print",& + description="Section of possible print options in FORCE_MIXING.",& + n_keywords=0, n_subsections=2, repeats=.FALSE., required=.TRUE.,& + error=error) + NULLIFY(print_key) + CALL cp_print_key_section_create(print_key,"SUBCELL",& + description="Activates the printing of the subcells used for the"//& + "generation of neighbor lists.", unit_str="angstrom",& + print_level=high_print_level,filename="__STD_OUT__",& + error=error) + CALL section_add_subsection(subsection,print_key,error=error) + CALL section_release(print_key,error=error) + + CALL cp_print_key_section_create(print_key,"NEIGHBOR_LISTS",& + description="Activates the printing of the neighbor lists used"//& + " for the hysteretic region calculations.", & + print_level=high_print_level,filename="",unit_str="angstrom",& + error=error) + CALL section_add_subsection(subsection,print_key,error=error) + CALL section_release(print_key,error=error) + + CALL section_add_subsection(section,subsection,error=error) + CALL section_release(subsection,error=error) + + + END SUBROUTINE create_qmmm_force_mixing_section + ! ***************************************************************************** !> \brief Input section to create QM kinds sections !> \param error variable to control error logging, stopping,... diff --git a/src/input_cp2k_restarts.F b/src/input_cp2k_restarts.F index f8378ed9dc..5e4385d957 100644 --- a/src/input_cp2k_restarts.F +++ b/src/input_cp2k_restarts.F @@ -112,7 +112,8 @@ MODULE input_cp2k_restarts PUBLIC :: write_restart,& update_subsys,& section_velocity_val_set,& - write_restart_header + write_restart_header,& + update_input CONTAINS diff --git a/src/integrator.F b/src/integrator.F index 6997101ae5..6df849c7da 100644 --- a/src/integrator.F +++ b/src/integrator.F @@ -79,6 +79,7 @@ MODULE integrator update_particle_set USE physcon, ONLY: bohr,& femtoseconds + USE qmmm_main, ONLY: qmmm_update_force_mixing_env USE qmmm_util, ONLY: apply_qmmm_walls_reflective USE reftraj_types, ONLY: reftraj_type USE reftraj_util, ONLY: compute_msd_reftraj @@ -255,6 +256,9 @@ CONTAINS DEALLOCATE (pos,STAT=stat) CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) + ![ADAPT] update input structure with new coordinates, make new labels + CALL qmmm_update_force_mixing_env(md_env, force_env, force_env%root_section, error=error) + ! Update forces CALL force_env_calc_energy_force(force_env,error=error) @@ -461,6 +465,9 @@ CONTAINS deallocate_vel=.FALSE. END IF + ![ADAPT] update input structure with new coordinates, make new labels + CALL qmmm_update_force_mixing_env(md_env, force_env, force_env%root_section, error=error) + CALL force_env_calc_energy_force(force_env,error=error) ! Metadynamics @@ -627,6 +634,9 @@ CONTAINS CALL update_dealloc_tmp ( tmp, particle_set, shell_particle_set, & core_particle_set, para_env, shell_adiabatic, pos=.TRUE., error=error ) + ![ADAPT] update input structure with new coordinates, make new labels + CALL qmmm_update_force_mixing_env(md_env, force_env, force_env%root_section, error=error) + CALL force_env_calc_energy_force(force_env,error=error) ! Metadynamics @@ -825,6 +835,9 @@ CONTAINS shell_particle_set, core_particle_set, globenv, tmp=tmp, check=.TRUE., error=error) END IF + ![ADAPT] update input structure with new coordinates, make new labels + CALL qmmm_update_force_mixing_env(md_env, force_env, force_env%root_section, error=error) + ! Update forces CALL force_env_calc_energy_force(force_env,error=error) @@ -1355,6 +1368,8 @@ CONTAINS ! Update forces CALL force_env_set_cell(force_env, cell, error) + ![ADAPT] update input structure with new coordinates, make new labels + CALL qmmm_update_force_mixing_env(md_env, force_env, force_env%root_section, error=error) CALL force_env_calc_energy_force(force_env,error=error) ! Metadynamics @@ -1576,6 +1591,9 @@ CPSourceFileRef,& CALL force_env_set_cell(force_env, cell, error) END IF + ![ADAPT] update input structure with new coordinates, make new labels + CALL qmmm_update_force_mixing_env(md_env, force_env, force_env%root_section, error=error) + ! Task to perform on the reference trajectory ! Compute energy and forces CALL force_env_calc_energy_force(force_env,eval_energy_forces=reftraj_env%info%eval_EF,error=error) @@ -1797,6 +1815,8 @@ CPSourceFileRef,& ! Update forces (and stress) CALL force_env_set_cell(force_env, cell, error) + ![ADAPT] update input structure with new coordinates, make new labels + CALL qmmm_update_force_mixing_env(md_env, force_env, force_env%root_section, error=error) CALL force_env_calc_energy_force(force_env, error=error) ! Metadynamics @@ -2059,6 +2079,8 @@ CPSourceFileRef,& ! Update forces CALL force_env_set_cell(force_env, cell, error) + ![ADAPT] update input structure with new coordinates, make new labels + CALL qmmm_update_force_mixing_env(md_env, force_env, force_env%root_section, error=error) CALL force_env_calc_energy_force(force_env, error=error) ! Metadynamics @@ -2357,6 +2379,8 @@ CPSourceFileRef,& ! Update forces CALL force_env_set_cell(force_env, cell, error) + ![ADAPT] update input structure with new coordinates, make new labels + CALL qmmm_update_force_mixing_env(md_env, force_env, force_env%root_section, error=error) CALL force_env_calc_energy_force(force_env,error=error) ! Metadynamics @@ -2579,6 +2603,9 @@ CPSourceFileRef,& particle_set_respa(iparticle)%r = particle_set(iparticle)%r END DO + ![ADAPT] update input structure with new coordinates, make new labels + ! Not checked, and probably won't work, with RESPA. + CALL qmmm_update_force_mixing_env(md_env, force_env, force_env%root_section, error=error) ! Update forces CALL force_env_calc_energy_force(force_env%sub_force_env(1)%force_env,error=error) @@ -2611,6 +2638,9 @@ CPSourceFileRef,& ! Multiple time step (second part) ! Compute forces for respa force_env + ![ADAPT] update input structure with new coordinates, make new labels + ! Not checked, and probably won't work, with RESPA. + CALL qmmm_update_force_mixing_env(md_env, force_env, force_env%root_section, error=error) CALL force_env_calc_energy_force(force_env,error=error) ! Metadynamics diff --git a/src/qmmm_force_mixing.F b/src/qmmm_force_mixing.F new file mode 100644 index 0000000000..ed992331e0 --- /dev/null +++ b/src/qmmm_force_mixing.F @@ -0,0 +1,735 @@ +! ***************************************************************************** +!> \brief Routines used for force-mixing QM/MM calculations +!> \par History +!> 2.2012 created [noam] +!> \author Noam Bernstein +! ***************************************************************************** +MODULE qmmm_force_mixing + USE atomic_kind_types, ONLY: atomic_kind_type,& + get_atomic_kind,& + set_atomic_kind + USE cell_types, ONLY: cell_type,& + pbc + USE cp_subsys_types, ONLY: cp_subsys_get,& + cp_subsys_type + USE f77_blas + USE fist_neighbor_list_types, ONLY: fist_neighbor_deallocate,& + fist_neighbor_type + USE fist_neighbor_lists, ONLY: build_fist_neighbor_lists + USE force_env_types, ONLY: force_env_get,& + force_env_type + USE input_section_types, ONLY: & + section_vals_add_values, section_vals_duplicate, section_vals_get, & + section_vals_get_subs_vals, section_vals_get_subs_vals3, & + section_vals_remove_values, section_vals_type, section_vals_val_get, & + section_vals_val_set + USE kinds, ONLY: default_string_length,& + dp + USE memory_utilities, ONLY: reallocate + USE mol_new_list_types, ONLY: mol_new_list_type + USE molecule_types_new, ONLY: molecule_type + USE particle_list_types, ONLY: particle_list_type + USE particle_types, ONLY: particle_type + USE qmmm_ff_fist, ONLY: qmmm_ff_precond_only_qm + USE qmmm_types, ONLY: force_mixing_label_QM_core,& + force_mixing_label_QM_core_list,& + force_mixing_label_QM_extended,& + force_mixing_label_buffer,& + force_mixing_label_none,& + force_mixing_label_termination +#include "cp_common_uses.h" + + IMPLICIT NONE + PRIVATE + + LOGICAL, PRIVATE :: debug_this_module=.FALSE. + CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qmmm_force_mixing' + + PUBLIC :: setup_force_mixing_qmmm_sections, update_force_mixing_labels + + +!*** +CONTAINS + +! ***************************************************************************** +!> \param error variable to control error logging, stopping,... +!> see module cp_error_handling +!> \par History +!> 02.2012 created [noam] +!> \author Noam Bernstein +! ***************************************************************************** + SUBROUTINE update_force_mixing_labels(force_env, subsys, qmmm_section, labels_changed, error) + TYPE(force_env_type), POINTER :: force_env + TYPE(cp_subsys_type), POINTER :: subsys + TYPE(section_vals_type), POINTER :: qmmm_section + LOGICAL, OPTIONAL :: labels_changed + TYPE(cp_error_type), INTENT(inout) :: error + + CHARACTER(len=*), PARAMETER :: routineN = 'update_force_mixing_labels', & + routineP = moduleN//':'//routineN + + INTEGER :: i_rep_section, i_rep_val, ip, im, & + n_new, n_rep_section, & + n_rep_val, natoms + INTEGER, ALLOCATABLE :: new_full_labels(:), & + orig_full_labels(:) + INTEGER, POINTER :: cur_indices(:), & + cur_labels(:), & + mm_index_entry(:), & + new_indices(:), new_labels(:) + LOGICAL :: explicit + REAL(dp), ALLOCATABLE :: nearest_dist(:) + REAL(dp), POINTER :: r_buf(:), r_core(:), r_qm(:) + TYPE(cell_type), POINTER :: cell + TYPE(fist_neighbor_type), POINTER :: nlist + TYPE(mol_new_list_type), POINTER :: molecules + TYPE(molecule_type), DIMENSION(:), & + POINTER :: molecule_set + TYPE(particle_list_type), POINTER :: particles + TYPE(particle_type), DIMENSION(:), & + POINTER :: particle_set + TYPE(section_vals_type), POINTER :: force_mixing_section, & + qm_kind_section, & + restart_section + integer :: QM_extended_seed_min_label_val + logical :: QM_extended_seed_is_core_list + integer, pointer :: breakable_bonds(:,:), broken_bonds(:) + integer, pointer :: molecule_of_atom(:) + integer :: max_n_qm + +if (debug_this_module) print *, "BOB starting update_force_mixing_labels" + ! get cur indices, labels + force_mixing_section => section_vals_get_subs_vals3(qmmm_section,"FORCE_MIXING",error=error) + CALL get_force_mixing_indices(force_mixing_section, cur_indices, cur_labels, error=error) +if (debug_this_module) print *, "BOB got cur_indices ",cur_indices(1:n_new) +if (debug_this_module) print *, "BOB got cur_labels ",cur_labels(1:n_new) + + ! read from input + ![NB] breakable bonds will come from here, too + NULLIFY(r_core, r_qm, r_buf) + CALL section_vals_val_get(force_mixing_section,"R_CORE",r_vals=r_core,error=error) + CALL section_vals_val_get(force_mixing_section,"R_QM",r_vals=r_qm,error=error) + CALL section_vals_val_get(force_mixing_section,"QM_EXTENDED_SEED_IS_ONLY_CORE_LIST",& + l_val=QM_extended_seed_is_core_list,error=error) + CALL section_vals_val_get(force_mixing_section,"R_BUF",r_vals=r_buf,error=error) + CALL section_vals_val_get(force_mixing_section,"MAX_N_QM",i_val=max_n_qm,error=error) + ![NB] need to read real list from input + ! should be 2xN_bb integer arrays, with (1,:) indices of inside atoms, and (2,:) indices of outside atoms + ! maybe also breakable_bond_types, with _atomic numbers_ of inside/outside atoms? + ! separate lists for core/buffer? + allocate(breakable_bonds(2,1)) + breakable_bonds(1:2,1) = 0 + + ! get particles, molecules + NULLIFY(particles,molecules) + CALL cp_subsys_get(subsys=subsys, particles=particles,molecules_new=molecules, error=error) + particle_set => particles%els + molecule_set => molecules%els + + natoms = SIZE(particle_set) + + ! initialize new indices, labels, and new_full_labels + NULLIFY(new_indices, new_labels) + CALL reallocate(new_indices,1,SIZE(cur_indices)) + CALL reallocate(new_labels,1,SIZE(cur_labels)) + new_indices = 0 + new_labels = force_mixing_label_none + ALLOCATE(new_full_labels(natoms)) + new_full_labels = force_mixing_label_none + + ! neighbor list for various hysteretic distance calls + NULLIFY(cell) + CALL force_env_get(force_env, cell=cell, error=error) + NULLIFY(nlist) + CALL make_neighbor_list(force_mixing_section, subsys, cell, MAX(r_core(2), r_qm(2), r_buf(2)), nlist, error) + + ! reset core_list from QM_KIND + NULLIFY(mm_index_entry) + qm_kind_section => section_vals_get_subs_vals3(qmmm_section,"QM_KIND",error=error) + CALL section_vals_get(qm_kind_section,n_repetition=n_rep_section,error=error) + n_new = 0 + DO i_rep_section=1,n_rep_section + CALL section_vals_val_get(qm_kind_section,"MM_INDEX",i_rep_section=i_rep_section,n_rep_val=n_rep_val,error=error) + DO i_rep_val=1,n_rep_val + CALL section_vals_val_get(qm_kind_section,"MM_INDEX",i_rep_section=i_rep_section,i_rep_val=i_rep_val, & + i_vals=mm_index_entry, error=error) + DO ip=1, SIZE(mm_index_entry) + call add_new_label(mm_index_entry(ip), force_mixing_label_QM_core_list, n_new, new_indices, new_labels, & + new_full_labels, max_n_qm) + END DO ! ip + END DO ! i_rep_val + END DO ! i_rep_section + +if (debug_this_module) print *, "BOB core_list new_indices ",new_indices(1:n_new) +if (debug_this_module) print *, "BOB core_list new_labels ",new_labels(1:n_new) + + ! allocate and initialize full atom set labels for hysteretic loops + ALLOCATE(nearest_dist(natoms)) + + ! orig_full_labels is full array (natoms) with orig labels + ALLOCATE(orig_full_labels(natoms)) + orig_full_labels = force_mixing_label_none + orig_full_labels(cur_indices(:)) = cur_labels(:) + + ! hysteretically set QM core from QM_core_list and radii, whole molecule + ![NB] need to replace all the whole molecule stuff with pad to breakable bonds. not quite done + ! (need intra molecule bond info, which isn't available for QM molecules yet) + + ! find molecule of atom + allocate(molecule_of_atom(natoms)) + do im=1, size(molecule_set) + molecule_of_atom(molecule_set(im)%first_atom:molecule_set(im)%last_atom) = im + end do + + ! add core using hysteretic selection(core_list, r_core) + unbreakable bonds + CALL add_layer_hysteretically(nlist, particle_set, cell, nearest_dist, & + orig_full_labels, new_full_labels, n_new, new_indices, new_labels, & + force_mixing_label_QM_core_list, force_mixing_label_QM_core_list, force_mixing_label_QM_core, r_core, & + max_n_qm, .TRUE., breakable_bonds, molecule_of_atom, molecule_set, broken_bonds) + ![NB] should actually pass this back for making link sections? + deallocate(broken_bonds) + +if (debug_this_module) print *, "BOB core new_indices ",new_indices(1:n_new) +if (debug_this_module) print *, "BOB core new_labels ",new_labels(1:n_new) + + ![NB] need more sophisticated QM extended, buffer rules + + ! add QM using hysteretic selection (core_list, r_qm) + unbreakable bonds +if (debug_this_module) print *, "BOB QM_extended_seed_is_core_list ",QM_extended_seed_is_core_list + if (QM_extended_seed_is_core_list) then + QM_extended_seed_min_label_val = force_mixing_label_QM_core_list + else ! QM region seed is all of core, not just core list + unbreakable bonds + QM_extended_seed_min_label_val = force_mixing_label_QM_core + endif + CALL add_layer_hysteretically(nlist, particle_set, cell, nearest_dist, & + orig_full_labels, new_full_labels, n_new, new_indices, new_labels, & + QM_extended_seed_min_label_val, force_mixing_label_QM_core_list, & + force_mixing_label_QM_extended, r_qm, & + max_n_qm, .TRUE., breakable_bonds, molecule_of_atom, molecule_set) + +if (debug_this_module) print *, "BOB extended new_indices ",new_indices(1:n_new) +if (debug_this_module) print *, "BOB extended new_labels ",new_labels(1:n_new) + + ! add buffer using hysteretic selection (>= QM extended, r_buf) + unbreakable bonds + CALL add_layer_hysteretically(nlist, particle_set, cell, nearest_dist, & + orig_full_labels, new_full_labels, n_new, new_indices, new_labels, & + force_mixing_label_QM_extended, force_mixing_label_QM_core_list, force_mixing_label_buffer, r_buf, & + max_n_qm, .TRUE., breakable_bonds, molecule_of_atom, molecule_set, broken_bonds) + ![NB] should actually pass this back for making link sections? + deallocate(broken_bonds) + +if (debug_this_module) print *, "BOB buffer new_indices ",new_indices(1:n_new) +if (debug_this_module) print *, "BOB buffer new_labels ",new_labels(1:n_new) + + deallocate(molecule_of_atom, breakable_bonds) + DEALLOCATE(nearest_dist) + + IF (PRESENT(labels_changed)) labels_changed = ANY (new_full_labels /= orig_full_labels) + + DEALLOCATE(orig_full_labels) + DEALLOCATE(new_full_labels) + + ! reduce new indices, labels to actually used size + CALL reallocate(new_indices,1,n_new) + CALL reallocate(new_labels,1,n_new) + + ! save info in input structure + restart_section => section_vals_get_subs_vals(qmmm_section, "FORCE_MIXING%RESTART_INFO", error=error) + CALL section_vals_get(restart_section,explicit=explicit,error=error) + IF (explicit) CALL section_vals_remove_values(restart_section,error) + CALL section_vals_val_set(restart_section,"INDICES",i_vals_ptr=new_indices,error=error) + CALL section_vals_val_set(restart_section,"LABELS",i_vals_ptr=new_labels,error=error) + + DEALLOCATE(cur_indices, cur_labels) + CALL fist_neighbor_deallocate(nlist,error) + + ![NB] perhap be controlled by some &PRINT section? + if (force_env%para_env%ionode) then + print *, "QMMM FORCE MIXING final count (not including links): "//& + " N_QM core ", count(new_labels >= force_mixing_label_QM_core), & + " N_QM buffered ", count(new_labels >= force_mixing_label_buffer) + endif + + END SUBROUTINE update_force_mixing_labels + + SUBROUTINE add_new_label(ip, label, n_new, new_indices, new_labels, new_full_labels, max_n_qm) + integer :: ip, label, n_new + integer, pointer :: new_indices(:), new_labels(:) + integer :: new_full_labels(:) + integer :: max_n_qm + + CHARACTER(len=*), PARAMETER :: routineN = 'update_force_mixing_labels', & + routineP = moduleN//':'//routineN + + n_new = n_new + 1 + CALL cp_assert(n_new <= max_n_qm, cp_failure_level,cp_assertion_failed,routineP,& + "add_new_label tried to add more atoms than allowed by &FORCE_MIXING&MAX_N_QM!"//& +CPSourceFileRef) + IF (n_new > SIZE(new_indices)) CALL reallocate(new_indices,1,n_new+9) + IF (n_new > SIZE(new_labels)) CALL reallocate(new_labels,1,n_new+9) + new_indices(n_new) = ip + new_labels(n_new) = label + new_full_labels(ip) = label + END SUBROUTINE add_new_label + + + SUBROUTINE add_layer_hysteretically(nlist, particle_set, cell, nearest_dist, & + orig_full_labels, new_full_labels, n_new, new_indices, new_labels, & + seed_min_label_val, seed_max_label_val, set_label_val, r_inout, max_n_qm, & + pad_to_breakable_bonds, breakable_bonds, molecule_of_atom, molecule_set, & + broken_bonds) + TYPE(fist_neighbor_type), POINTER :: nlist + TYPE(particle_type), DIMENSION(:), & + POINTER :: particle_set + TYPE(cell_type), POINTER :: cell + REAL(dp) :: nearest_dist(:) + INTEGER :: orig_full_labels(:), new_full_labels(:) + integer :: n_new + INTEGER, POINTER :: new_indices(:), new_labels(:) + INTEGER :: seed_min_label_val, & + seed_max_label_val, & + set_label_val + REAL(dp) :: r_inout(2) + integer :: max_n_qm + LOGICAL :: pad_to_breakable_bonds + integer, pointer, optional :: molecule_of_atom(:), breakable_bonds(:,:), broken_bonds(:) + integer :: bond_atom_a, bond_atom_b, bond_atom_inside, bond_atom_outside + TYPE(molecule_type), DIMENSION(:), & + OPTIONAL, POINTER :: molecule_set + + integer :: natoms + + CHARACTER(len=*), PARAMETER :: routineN = 'add_layer_hysteretically', & + routineP = moduleN//':'//routineN + + INTEGER :: i_ind, im, ip, ib, ibb, ipair, & + ipairkind, j_ind + LOGICAL :: molec_in_inner, molec_in_outer + LOGICAL :: i_in_new_seed, j_in_new_seed + LOGICAL :: i_outside_new_seed, j_outside_new_seed + LOGICAL :: breakable + REAL(dp) :: r_ij(3), r_ij_mag + integer :: orig_n_new + logical :: any_additions + integer :: n_broken_bonds + +if (debug_this_module) print *,"BOB adding hysteretically seed ",seed_min_label_val,seed_max_label_val,& +" set ",set_label_val," r ",r_inout + ! calculate nearest dist from each atom outside of new seed to nearest atom inside of new seed + nearest_dist=HUGE(1.0_dp) + ! loop over pairs of all kinds in random order + DO ipairkind=1,SIZE(nlist%neighbor_kind_pairs) + DO ipair=1,nlist%neighbor_kind_pairs(ipairkind)%npairs + + i_ind = nlist%neighbor_kind_pairs(ipairkind)%list(1,ipair) + j_ind = nlist%neighbor_kind_pairs(ipairkind)%list(2,ipair) + + i_in_new_seed = (new_full_labels(i_ind) >= seed_min_label_val .AND. new_full_labels(i_ind) <= seed_max_label_val) + i_outside_new_seed = (new_full_labels(i_ind) < seed_min_label_val) + j_in_new_seed = (new_full_labels(j_ind) >= seed_min_label_val .AND. new_full_labels(j_ind) <= seed_max_label_val) + j_outside_new_seed = (new_full_labels(j_ind) < seed_min_label_val) + + IF ((i_in_new_seed .and. j_outside_new_seed) .or. (j_in_new_seed .and. i_outside_new_seed)) THEN + r_ij = pbc(particle_set(i_ind)%r - particle_set(j_ind)%r, cell) + r_ij_mag = SQRT(SUM(r_ij**2)) + IF (i_in_new_seed .and. j_outside_new_seed .and. (r_ij_mag < nearest_dist(j_ind))) then + nearest_dist(j_ind) = r_ij_mag + ENDIF + IF (j_in_new_seed .and. i_outside_new_seed .and. (r_ij_mag < nearest_dist(i_ind))) then + nearest_dist(i_ind) = r_ij_mag + ENDIF + ENDIF + + END DO + END DO +if (debug_this_module) print *,"BOB nearest_dist ", nearest_dist + +![NB] this is whole molecule. Should be replaced with labeling of individual atoms + +! pad_to_breakable_bonds (below), but QM molecule bond information isn't available yet +DO im=1, SIZE(molecule_set) +if (debug_this_module) print *, "BOB check molecule ",im," atoms ",molecule_set(im)%first_atom,molecule_set(im)%last_atom + molec_in_inner = any(nearest_dist(molecule_set(im)%first_atom:molecule_set(im)%last_atom) <= r_inout(1)) + molec_in_outer = any(nearest_dist(molecule_set(im)%first_atom:molecule_set(im)%last_atom) <= r_inout(2)) +if (debug_this_module) print *, "BOB got molec_in_inner ", molec_in_inner, " outer ", molec_in_outer + if (molec_in_inner) then + do ip=molecule_set(im)%first_atom, molecule_set(im)%last_atom + ! labels are being rebuild from scratch, so never overwrite new label that's higher level + if (new_full_labels(ip) < set_label_val) & + call add_new_label(ip, set_label_val, n_new, new_indices, new_labels, new_full_labels, max_n_qm) + end do + else if (molec_in_outer) then + if (any(orig_full_labels(molecule_set(im)%first_atom:molecule_set(im)%last_atom) >= set_label_val)) then + do ip=molecule_set(im)%first_atom, molecule_set(im)%last_atom + ! labels are being rebuild from scratch, so never overwrite new label that's higher level + if (new_full_labels(ip) < set_label_val) & + call add_new_label(ip, set_label_val, n_new, new_indices, new_labels, new_full_labels, max_n_qm) + end do + endif + endif +end do +if (present(broken_bonds)) CALL reallocate(broken_bonds,1,0) + +!!![NB] enable this (and remove previous bit) once intra-QM-molecule bond information +!!! is available +!!! ! label new layer +!!! ! loop over every atom +!!! DO ip=1, SIZE(new_full_labels) +!!! ! labels are being rebuild from scratch, so never overwrite new label that's higher level +!!! if (new_full_labels(ip) > set_label_val) cycle +!!! +!!!if (debug_this_module) print *,"BOB atom ",ip," new_full_label ",new_full_labels(ip), & +!!!" is not in seed, check for for adding" +!!! ![NB] is this logic right? maybe it is +!!! IF ( nearest_dist(ip) <= r_inout(1) .OR. & +!!! (nearest_dist(ip) <= r_inout(2) .AND. orig_full_labels(ip) >= set_label_val) ) THEN +!!! ! atom in inner radius +!!! ! or +!!! ! atom in outer radius and was previously set to set_label_val or higher +!!!if (debug_this_module) print *,"BOB add label for atom ", ip +!!! call add_new_label(ip, set_label_val, n_new, new_indices, new_labels, new_full_labels, max_n_qm) +!!! ENDIF +!!! END DO ! ip +!!! +!!! if (present(broken_bonds)) n_broken_bonds = 0 +!!! if (pad_to_breakable_bonds) then +!!!if (debug_this_module) print *,"BOB pad to breakable bonds" +!!! CALL cp_assert(PRESENT(molecule_of_atom),cp_failure_level,cp_assertion_failed,routineP,& +!!! "pad_to_breakable_bonds is true, but molecule_of_atom not present!"//& +!!!CPSourceFileRef) +!!! CALL cp_assert(PRESENT(molecule_set),cp_failure_level,cp_assertion_failed,routineP,& +!!! "pad_to_breakable_bonds is true, but molecule_set not present!"//& +!!!CPSourceFileRef) +!!! CALL cp_assert(PRESENT(breakable_bonds),cp_failure_level,cp_assertion_failed,routineP,& +!!! "whole_molecule is true, but breakable_bonds not present!"//& +!!!CPSourceFileRef) +!!! +!!! orig_n_new = n_new +!!! ! loop over all recently marked atoms +!!! do ip=1, orig_n_new +!!! if (new_labels(ip) /= set_label_val) cycle +!!! +!!!if (debug_this_module) then +!!! print *, "BOB molecule set" +!!! do im=1, size(molecule_set) +!!! print *, "BOB molecule im ",im, " #bonds ",size(molecule_set(im)%molecule_kind%bond_list), " name ", & +!!! trim(molecule_set(im)%molecule_kind%name), molecule_set(im)%molecule_kind%kind_number +!!! do ib=1, size(molecule_set(im)%molecule_kind%bond_list) +!!! print *, "BOB",molecule_set(im)%molecule_kind%bond_list(ib)%a + molecule_set(im)%first_atom - 1, & +!!! molecule_set(im)%molecule_kind%bond_list(ib)%b + molecule_set(im)%first_atom - 1 +!!! end do +!!! end do +!!!endif +!!! +!!!if (debug_this_module) print *,"BOB check atom local index ",ip," full index ", new_indices(ip) +!!! ! atom has just been set to set_label_val, check the rest of its molecule +!!! im = molecule_of_atom(new_indices(ip)) +!!! +!!!if (debug_this_module) print *,"in molecule ",im +!!! any_additions = .true. +!!! do while (any_additions) +!!!if (debug_this_module) print *,"start loop checking for additions" +!!! any_additions = .false. +!!! ! loop over bonds in this molecule in random order +!!! do ib=1, size(molecule_set(im)%molecule_kind%bond_list) +!!! bond_atom_a = molecule_set(im)%molecule_kind%bond_list(ib)%a + molecule_set(im)%first_atom - 1 +!!! bond_atom_b = molecule_set(im)%molecule_kind%bond_list(ib)%b + molecule_set(im)%first_atom - 1 +!!!if (debug_this_module) print *,"bond ", ib, " between a,b ", bond_atom_a,bond_atom_b +!!! +!!! if (new_full_labels(bond_atom_a) == set_label_val .and. new_full_labels(bond_atom_b) < set_label_val) then +!!! ! a is inside, b is outside +!!! bond_atom_inside = bond_atom_a +!!! bond_atom_outside = bond_atom_b +!!!if (debug_this_module) print *," a is inside, b is outside" +!!! else if (new_full_labels(bond_atom_b) == set_label_val .and. new_full_labels(bond_atom_a) < set_label_val) then +!!! ! b is inside, a is outside +!!! bond_atom_inside = bond_atom_b +!!! bond_atom_outside = bond_atom_a +!!!if (debug_this_module) print *," b is inside, a is outside" +!!! else +!!! ! inside-inside or outside-outside +!!! cycle +!!! endif +!!! +!!! ! check for breakable bonds +!!! breakable = .false. +!!! do ibb=1,size(breakable_bonds,2) +!!! if ((breakable_bonds(1,ibb) == bond_atom_inside) .and. breakable_bonds(2,ibb) == bond_atom_outside) then +!!! ! this bond is OK to break. +!!! breakable = .true. +!!! if (set_label_val == force_mixing_label_buffer) then ! set termination mark on outside atom b +!!! call add_new_label(bond_atom_outside, force_mixing_label_termination, n_new, new_indices, new_labels, & +!!! new_full_labels, max_n_qm) +!!! endif +!!! if (present(broken_bonds)) then ! add a new one +!!! n_broken_bonds = n_broken_bonds + 1 +!!! CALL reallocate(broken_bonds,1,n_broken_bonds) +!!! broken_bonds(n_broken_bonds) = ibb +!!! endif +!!! endif +!!! end do ! ibb +!!!if (debug_this_module) print *," breakable ", breakable +!!! +!!! if (.not. breakable) then ! add correct atom +!!! call add_new_label(bond_atom_outside, set_label_val, n_new, new_indices, new_labels, new_full_labels, max_n_qm) +!!! any_additions = .true. +!!!if (debug_this_module) print *, " not breakable, adding outside atom ", bond_atom_outside +!!! endif +!!! +!!! end do ! ib +!!! end do ! any_additions +!!! +!!! end do ! ip +!!! +!!! if (present(broken_bonds)) CALL reallocate(broken_bonds,1,n_broken_bonds) +!!! endif ! pad_to_breakable_bonds + + END SUBROUTINE add_layer_hysteretically + + SUBROUTINE make_neighbor_list(force_mixing_section, subsys, cell, r_max, nlist, error) + TYPE(section_vals_type), POINTER :: force_mixing_section + TYPE(cp_subsys_type), POINTER :: subsys + TYPE(cell_type), POINTER :: cell + REAL(dp) :: r_max + TYPE(fist_neighbor_type), POINTER :: nlist + TYPE(cp_error_type), INTENT(inout) :: error + + CHARACTER(len=*), PARAMETER :: routineN = 'make_neighbor_list', & + routineP = moduleN//':'//routineN + + CHARACTER(LEN=default_string_length) :: kind_name + CHARACTER(LEN=default_string_length), & + POINTER :: kind_name_a(:) + INTEGER :: ik + LOGICAL :: skip_kind + REAL(dp), ALLOCATABLE :: r_max_a(:,:), r_minsq_a(:,:) + TYPE(atomic_kind_type), POINTER :: atomic_kind + + ALLOCATE(r_max_a(SIZE(subsys%atomic_kinds%els),SIZE(subsys%atomic_kinds%els))) + ALLOCATE(r_minsq_a(SIZE(subsys%atomic_kinds%els),SIZE(subsys%atomic_kinds%els))) + r_max_a = r_max + r_minsq_a = EPSILON(1.0_dp) + + ! save kind names + ALLOCATE(kind_name_a(SIZE(subsys%atomic_kinds%els))) + DO ik=1, SIZE(subsys%atomic_kinds%els) + atomic_kind => subsys%atomic_kinds%els(ik) + CALL get_atomic_kind(atomic_kind=atomic_kind,name=kind_name) + kind_name_a(ik) = kind_name + END DO + + ! overwrite kind names so that none are QM, and so excluding QM-QM interactions + ! (which is not what we want) will not happen + DO ik=1, SIZE(subsys%atomic_kinds%els) + atomic_kind => subsys%atomic_kinds%els(ik) + CALL get_atomic_kind(atomic_kind=atomic_kind,name=kind_name) + ! when atom is QM atom, kind_name is replaced with original + ! mm kind name, and return status is logical .TRUE. + skip_kind=qmmm_ff_precond_only_qm(kind_name) + CALL set_atomic_kind(atomic_kind=atomic_kind,name=kind_name) + END DO + + NULLIFY(nlist) + CALL build_fist_neighbor_lists(subsys%atomic_kinds%els, subsys%particles%els, & + cell=cell, r_max=r_max_a, r_minsq=r_minsq_a, & + ei_scale14=1.0_dp, vdw_scale14=1.0_dp, nonbonded=nlist, & + para_env=subsys%para_env, build_from_scratch=.TRUE., geo_check=.FALSE., & + mm_section=force_mixing_section, error=error) + + DEALLOCATE(r_max_a, r_minsq_a) + + ! restore kind names + DO ik=1, SIZE(subsys%atomic_kinds%els) + CALL set_atomic_kind(atomic_kind=atomic_kind,name=kind_name_a(ik)) + END DO + DEALLOCATE(kind_name_a) + + END SUBROUTINE make_neighbor_list + +! ***************************************************************************** +!> \param error variable to control error logging, stopping,... +!> see module cp_error_handling +!> \par History +!> 02.2012 created [noam] +!> \author Noam Bernstein +! ***************************************************************************** + SUBROUTINE setup_force_mixing_qmmm_sections(subsys, qmmm_section, qmmm_core_section, qmmm_extended_section, error) + TYPE(cp_subsys_type), POINTER :: subsys + TYPE(section_vals_type), POINTER :: qmmm_section, & + qmmm_core_section, & + qmmm_extended_section + TYPE(cp_error_type), INTENT(inout) :: error + + CHARACTER(len=*), PARAMETER :: & + routineN = 'setup_force_mixing_qmmm_sections', & + routineP = moduleN//':'//routineN + + CHARACTER(len=DEFAULT_STRING_LENGTH), & + POINTER :: elem_mapping(:,:), & + elem_mapping_entry(:) + INTEGER :: i_rep_section_core, i_rep_section_extended, i_rep_val_core, & + i_rep_val_extended, ielem, ip, n_elements + INTEGER, POINTER :: cur_indices(:), cur_labels(:) + LOGICAL :: new_element_core, & + new_element_extended + TYPE(particle_type), DIMENSION(:), & + POINTER :: particles + TYPE(section_vals_type), POINTER :: force_mixing_section, & + link_section, qm_kind_section + +! +! create new qmmm sections for core and extended + + NULLIFY(qmmm_core_section, qmmm_extended_section) + CALL section_vals_duplicate(qmmm_section, qmmm_core_section, error=error) + CALL section_vals_duplicate(qmmm_section, qmmm_extended_section, error=error) + + ! remove LINKs (specified by user for core) from extended + link_section => section_vals_get_subs_vals(qmmm_extended_section, "LINK", can_return_null=.TRUE., error=error) + IF (ASSOCIATED(link_section)) THEN + CALL section_vals_remove_values(link_section,error) + END IF + + force_mixing_section => section_vals_get_subs_vals(qmmm_section,"FORCE_MIXING",error=error) + + ! get QM_KIND_ELEMENT_MAPPING + CALL section_vals_val_get(force_mixing_section,"QM_KIND_ELEMENT_MAPPING",n_rep_val=n_elements,error=error) + ALLOCATE(elem_mapping(2,n_elements)) + DO ielem=1,n_elements + CALL section_vals_val_get(force_mixing_section,"QM_KIND_ELEMENT_MAPPING",i_rep_val=ielem,c_vals=elem_mapping_entry,& + error=error) + elem_mapping(1:2,ielem) = elem_mapping_entry(1:2) + END DO + + ! get CUR_INDICES, CUR_LABELS + CALL get_force_mixing_indices(force_mixing_section, cur_indices, cur_labels, error=error) + CALL cp_assert(SIZE(cur_indices) > 0,cp_failure_level,cp_assertion_failed,routineP,& + "cur_indices is empty, found no QM atoms"//& +CPSourceFileRef) + + ! loop through elements and atoms, and set up new QM_KIND sections + particles => subsys%particles%els + + qm_kind_section => section_vals_get_subs_vals3(qmmm_section,"QM_KIND",error=error) + CALL section_vals_get(qm_kind_section,n_repetition=i_rep_section_core,error=error) + CALL cp_assert(i_rep_section_core > 0,cp_failure_level,cp_assertion_failed,routineP,& + "Force-mixing QM didn't find any QM_KIND sections, so no core specified!"//& +CPSourceFileRef) + i_rep_section_extended = i_rep_section_core + DO ielem=1,n_elements + new_element_core = .TRUE. + new_element_extended = .TRUE. + DO ip=1, SIZE(cur_indices) ! particles with label + + ! extended + ! if current particle is some sort of QM atom, and not in core list (those the user + ! gave explicit QM_KIND sections for, and not a termination atom, need to make a + ! QM_KIND section for it + IF (cur_labels(ip) > force_mixing_label_none .AND. cur_labels(ip) < force_mixing_label_QM_core_list .AND. & + cur_labels(ip) /= force_mixing_label_termination) THEN + IF (TRIM(particles(cur_indices(ip))%atomic_kind%element_symbol) == TRIM(elem_mapping(1,ielem))) THEN + qm_kind_section => section_vals_get_subs_vals3(qmmm_extended_section,"QM_KIND",error=error) + IF (new_element_extended) THEN ! add new QM_KIND section for this element + i_rep_section_extended = i_rep_section_extended + 1 + CALL section_vals_add_values(qm_kind_section,error=error) + CALL section_vals_val_set(qm_kind_section,"_SECTION_PARAMETERS_",i_rep_section=i_rep_section_extended, & + c_val=elem_mapping(2,ielem),error=error) + i_rep_val_extended = 0 + new_element_extended = .FALSE. + ENDIF + i_rep_val_extended = i_rep_val_extended + 1 + CALL section_vals_val_set(qm_kind_section,"MM_INDEX",i_rep_section=i_rep_section_extended, & + i_rep_val=i_rep_val_extended, i_val=cur_indices(ip),error=error) + ENDIF ! element name matches current element + ENDIF ! is a non-termination QM atom + + ! core + ! if current particle is a core QM atom, and not in core list (those the user + ! gave explicit QM_KIND sections for, need to make a QM_KIND section for it + IF (cur_labels(ip) == force_mixing_label_QM_core) THEN + IF (TRIM(particles(cur_indices(ip))%atomic_kind%element_symbol) == TRIM(elem_mapping(1,ielem))) THEN + qm_kind_section => section_vals_get_subs_vals3(qmmm_core_section,"QM_KIND",error=error) + IF (new_element_core) THEN ! add new QM_KIND section for this element + i_rep_section_core = i_rep_section_core + 1 + CALL section_vals_add_values(qm_kind_section,error=error) + CALL section_vals_val_set(qm_kind_section,"_SECTION_PARAMETERS_",i_rep_section=i_rep_section_core, & + c_val=elem_mapping(2,ielem),error=error) + i_rep_val_core = 0 + new_element_core = .FALSE. + ENDIF + i_rep_val_core = i_rep_val_core + 1 + CALL section_vals_val_set(qm_kind_section,"MM_INDEX",i_rep_section=i_rep_section_core, & + i_rep_val=i_rep_val_core, i_val=cur_indices(ip),error=error) + ENDIF ! element name matches current element + ENDIF ! is a non-termination QM atom + + END DO ! atom index ip + END DO ! element index ielem + + ![NB] check + DEALLOCATE(elem_mapping, cur_indices, cur_labels) + + END SUBROUTINE setup_force_mixing_qmmm_sections + + SUBROUTINE get_force_mixing_indices(force_mixing_section, indices, labels, error) + TYPE(section_vals_type), POINTER :: force_mixing_section + INTEGER, POINTER :: indices(:), labels(:) + TYPE(cp_error_type), INTENT(inout) :: error + + CHARACTER(len=*), PARAMETER :: routineN = 'get_force_mixing_indices', & + routineP = moduleN//':'//routineN + + INTEGER :: i_rep_val, n_indices, & + n_labels, n_reps + INTEGER, POINTER :: indices_entry(:), & + labels_entry(:) + LOGICAL :: explicit + TYPE(section_vals_type), POINTER :: restart_section + + NULLIFY(indices, labels) + restart_section => section_vals_get_subs_vals(force_mixing_section, "RESTART_INFO", error=error) + CALL section_vals_get(restart_section,explicit=explicit,error=error) + IF (.NOT. explicit) THEN ! no old indices, labels, return empty arrays + ALLOCATE(indices(0)) + ALLOCATE(labels(0)) + RETURN + ENDIF + + ![NB] maybe switch to reallocatable array + CALL section_vals_val_get(restart_section,"INDICES",n_rep_val=n_reps,error=error) + n_indices=0 + DO i_rep_val = 1,n_reps + CALL section_vals_val_get(restart_section,"INDICES",& + i_rep_val=i_rep_val,i_vals=indices_entry,error=error) + n_indices = n_indices + SIZE(indices_entry) + END DO + ALLOCATE(indices(n_indices)) + n_indices=0 + DO i_rep_val = 1,n_reps + CALL section_vals_val_get(restart_section,"INDICES",& + i_rep_val=i_rep_val,i_vals=indices_entry,error=error) + indices(n_indices+1:n_indices+SIZE(indices_entry)) = indices_entry + n_indices = n_indices + SIZE(indices_entry) + END DO + + CALL section_vals_val_get(restart_section,"LABELS",n_rep_val=n_reps,error=error) + n_labels=0 + DO i_rep_val = 1,n_reps + CALL section_vals_val_get(restart_section,"LABELS",& + i_rep_val=i_rep_val,i_vals=labels_entry,error=error) + n_labels = n_labels + SIZE(labels_entry) + END DO + ALLOCATE(labels(n_labels)) + n_labels=0 + DO i_rep_val = 1,n_reps + CALL section_vals_val_get(restart_section,"LABELS",& + i_rep_val=i_rep_val,i_vals=labels_entry,error=error) + labels(n_labels+1:n_labels+SIZE(labels_entry)) = labels_entry + n_labels = n_labels + SIZE(labels_entry) + END DO + + CALL cp_assert(n_indices == n_labels, cp_failure_level,cp_assertion_failed,routineP,& + "got unequal numbers of force_mixing indices and labels!"//& +CPSourceFileRef) + END SUBROUTINE get_force_mixing_indices + +END MODULE qmmm_force_mixing diff --git a/src/qmmm_main.F b/src/qmmm_main.F index 83393d1983..d338fab0ad 100644 --- a/src/qmmm_main.F +++ b/src/qmmm_main.F @@ -23,6 +23,7 @@ MODULE qmmm_main USE cp_para_types, ONLY: cp_para_env_type USE cp_subsys_methods, ONLY: create_small_subsys USE cp_subsys_types, ONLY: cp_subsys_release,& + cp_subsys_retain,& cp_subsys_type USE f77_blas USE fist_environment_types, ONLY: fist_env_get,& @@ -32,6 +33,7 @@ MODULE qmmm_main USE force_env_methods, ONLY: force_env_create USE force_env_types, ONLY: force_env_get,& force_env_p_type,& + force_env_release,& force_env_type USE global_types, ONLY: global_environment_type USE header, ONLY: qmmm_header @@ -39,13 +41,18 @@ MODULE qmmm_main do_qmmm,& do_qmmm_none,& do_qs + USE input_cp2k_restarts, ONLY: update_input USE input_section_types, ONLY: section_vals_get,& section_vals_get_subs_vals,& + section_vals_release,& section_vals_type,& section_vals_val_get,& section_vals_val_set USE kinds, ONLY: default_string_length,& dp + USE md_environment_types, ONLY: md_environment_type + USE qmmm_force_mixing, ONLY: setup_force_mixing_qmmm_sections,& + update_force_mixing_labels USE qmmm_init, ONLY: & assign_mm_charges_and_radius, move_or_add_atoms, & print_image_charge_info, print_qmmm_charges, print_qmmm_links, & @@ -55,15 +62,19 @@ MODULE qmmm_main USE qmmm_links_methods, ONLY: qmmm_link_Imomm_coord USE qmmm_pw_grid, ONLY: qmmm_pw_grid_init USE qmmm_types, ONLY: & - add_set_release, add_set_type, qmmm_env_mm_create, & + add_set_release, add_set_type, fist_subsys, force_mixing_core_subsys, & + force_mixing_extended_subsys, primary_subsys, qmmm_env_mm_create, & qmmm_env_mm_release, qmmm_env_mm_type, qmmm_env_qm_create, & - qmmm_env_qm_release, qmmm_env_qm_type, qmmm_links_type + qmmm_env_qm_release, qmmm_env_qm_retain, qmmm_env_qm_type, & + qmmm_links_type, qs_subsys USE qs_environment_types, ONLY: qs_environment_type,& set_qs_env USE qs_main, ONLY: quickstep_create_force_env USE termination, ONLY: stop_program USE timings, ONLY: timeset,& timestop + USE virial_types, ONLY: virial_release,& + virial_retain #include "cp_common_uses.h" IMPLICIT NONE @@ -72,10 +83,182 @@ MODULE qmmm_main LOGICAL, PRIVATE, PARAMETER :: debug_this_module=.TRUE. CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qmmm_main' - PUBLIC :: qmmm_create_force_env + PUBLIC :: qmmm_create_force_env, qmmm_update_force_mixing_env + + !*** CONTAINS +! ***************************************************************************** +!> \param error variable to control error logging, stopping,... +!> see module cp_error_handling +!> \par History +!> 02.2012 created [noam] +!> \author Noam Bernstein +! ***************************************************************************** + SUBROUTINE qmmm_create_force_env(force_env, root_section, para_env, globenv,& + force_env_section, subsys_section, use_motion_section, dummy_force_env, error) + TYPE(force_env_type), POINTER :: force_env + TYPE(section_vals_type), POINTER :: root_section + TYPE(cp_para_env_type), POINTER :: para_env + TYPE(global_environment_type), POINTER :: globenv + TYPE(section_vals_type), POINTER :: force_env_section, & + subsys_section + LOGICAL, INTENT(IN) :: use_motion_section + TYPE(force_env_type), OPTIONAL, POINTER :: dummy_force_env + TYPE(cp_error_type), INTENT(inout) :: error + + CHARACTER(len=*), PARAMETER :: routineN = 'qmmm_create_force_env', & + routineP = moduleN//':'//routineN + + INTEGER :: stat + LOGICAL :: failure, force_mixing_active + TYPE(cp_subsys_type), POINTER :: subsys + TYPE(force_env_p_type), DIMENSION(:), & + POINTER :: sub_force_env + TYPE(force_env_type), POINTER :: use_dummy_force_env + TYPE(qmmm_env_qm_type), POINTER :: qmmm_env_qm + TYPE(section_vals_type), POINTER :: qmmm_core_section, & + qmmm_extended_section, & + qmmm_force_mixing, & + qmmm_section + + qmmm_section => section_vals_get_subs_vals(force_env_section,"QMMM",error=error) + qmmm_force_mixing => section_vals_get_subs_vals(qmmm_section,"FORCE_MIXING",error=error) + CALL section_vals_get(qmmm_force_mixing,explicit=force_mixing_active,error=error) + + IF (force_mixing_active) THEN + IF (PRESENT(dummy_force_env)) THEN + use_dummy_force_env => force_env + ELSE + CALL qmmm_create_force_env_low(use_dummy_force_env, & + root_section, para_env, globenv,& + force_env_section, qmmm_section, subsys_section, use_motion_section, error) + ENDIF + CALL force_env_get(use_dummy_force_env, subsys=subsys, error=error) + + ![NB] + ! check that nothing accesses via root or force_env sections + ALLOCATE(sub_force_env(2),stat=stat) + CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) + + CALL update_force_mixing_labels(use_dummy_force_env, subsys, qmmm_section, error=error) + + ! using CUR_INDICES and CUR_LABELS, create appropriate QM_KIND sections for two QM/MM calculations + CALL setup_force_mixing_qmmm_sections(subsys, qmmm_section, qmmm_core_section, qmmm_extended_section,error=error) + + NULLIFY(sub_force_env(force_mixing_core_subsys)%force_env) + CALL qmmm_create_force_env_low(sub_force_env(force_mixing_core_subsys)%force_env, & + root_section, para_env, globenv,& + force_env_section, qmmm_core_section, subsys_section, use_motion_section, error) + NULLIFY(sub_force_env(force_mixing_extended_subsys)%force_env) + CALL qmmm_create_force_env_low(sub_force_env(force_mixing_extended_subsys)%force_env, & + root_section, para_env, globenv,& + force_env_section, qmmm_extended_section, subsys_section, use_motion_section, error) + + CALL section_vals_release(qmmm_core_section, error=error) + CALL section_vals_release(qmmm_extended_section, error=error) + + IF (.NOT. PRESENT(dummy_force_env)) THEN ! allocated it above, now release + CALL force_env_release(use_dummy_force_env, error=error) + ENDIF + ELSE + ALLOCATE(sub_force_env(1),stat=stat) + CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) + CALL qmmm_create_force_env_low(sub_force_env(1)%force_env, & + root_section, para_env, globenv,& + force_env_section, qmmm_section, subsys_section, use_motion_section, error) + ENDIF + + CALL force_env_get(sub_force_env(primary_subsys)%force_env,qmmm_env=qmmm_env_qm,error=error) + CALL force_env_create(force_env, root_section, sub_force_env=sub_force_env,& + qmmm_env=qmmm_env_qm, & + para_env=para_env, globenv=globenv, force_env_section= force_env_section,& + error=error) + + END SUBROUTINE qmmm_create_force_env + + + SUBROUTINE qmmm_update_force_mixing_env(md_env, force_env, root_section, error) + TYPE(md_environment_type), POINTER :: md_env + TYPE(force_env_type), POINTER :: force_env + TYPE(section_vals_type), POINTER :: root_section + TYPE(cp_error_type), INTENT(inout) :: error + + CHARACTER(len=*), PARAMETER :: routineN = 'qmmm_update_force_mixing_env', & + routineP = moduleN//':'//routineN + + INTEGER :: stat + LOGICAL :: failure, force_mixing_active, & + labels_changed + TYPE(cp_subsys_type), POINTER :: subsys + TYPE(force_env_p_type), DIMENSION(:), & + POINTER :: sub_force_env + TYPE(section_vals_type), POINTER :: qmmm_core_section, & + qmmm_extended_Section, qmmm_force_mixing, qmmm_section, subsys_section + + qmmm_section => section_vals_get_subs_vals(force_env%force_env_section,"QMMM",error=error) + qmmm_force_mixing => section_vals_get_subs_vals(qmmm_section,"FORCE_MIXING",error=error) + CALL section_vals_get(qmmm_force_mixing,explicit=force_mixing_active,error=error) + + IF (force_mixing_active) THEN + + CALL force_env_get(force_env, subsys=subsys, error=error) + CALL update_force_mixing_labels(force_env, subsys, qmmm_section, labels_changed=labels_changed, error=error) + IF (.NOT. labels_changed) THEN + RETURN + ENDIF + CALL cp_assert(.FALSE.,cp_warning_level,cp_assertion_failed,routineP,& + "Adaptive force-mixing labels changed, rebuilding QM/MM calculations! "//& +CPSourceFileRef,& + only_ionode=.TRUE.) + + CALL update_input(md_env,force_env,root_section,error=error) + + ALLOCATE(sub_force_env(2),stat=stat) + CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) + + ! using CUR_INDICES and CUR_LABELS, create appropriate QM_KIND sections for two QM/MM calculations + CALL setup_force_mixing_qmmm_sections(subsys, qmmm_section, qmmm_core_section, qmmm_extended_section,error=error) + + subsys_section => section_vals_get_subs_vals(force_env%force_env_section,"SUBSYS",error=error) + NULLIFY(sub_force_env(force_mixing_core_subsys)%force_env) + ![ADAPT] no sure about use_motion_section + CALL qmmm_create_force_env_low(sub_force_env(force_mixing_core_subsys)%force_env, & + force_env%root_section, force_env%para_env, force_env%globenv,& + force_env%force_env_section, qmmm_core_section, subsys_section, use_motion_section=.TRUE., error=error) + NULLIFY(sub_force_env(force_mixing_extended_subsys)%force_env) + CALL qmmm_create_force_env_low(sub_force_env(force_mixing_extended_subsys)%force_env, & + force_env%root_section, force_env%para_env, force_env%globenv,& + force_env%force_env_section, qmmm_extended_section, subsys_section, use_motion_section=.TRUE., error=error) + + CALL section_vals_release(qmmm_core_section, error=error) + CALL section_vals_release(qmmm_extended_section, error=error) + + ! do something with sub_force_envs + + ! duplicate releases in force_env_release + CALL force_env_release(force_env%sub_force_env(force_mixing_core_subsys)%force_env,error=error) + CALL force_env_release(force_env%sub_force_env(force_mixing_extended_subsys)%force_env,error=error) + CALL qmmm_env_qm_release(force_env%qmmm_env,error=error) + CALL virial_release(force_env%virial,error=error) + CALL cp_subsys_release(force_env%subsys,error=error) + + DEALLOCATE(force_env%sub_force_env) + force_env%sub_force_env => sub_force_env + + ! duplicate things that are pointed to in force_env_create when doing qmmm + force_env%virial => sub_force_env(primary_subsys)%force_env%virial + CALL virial_retain(force_env%virial, error=error) + force_env%qmmm_env => sub_force_env(primary_subsys)%force_env%qmmm_env + CALL qmmm_env_qm_retain(force_env%qmmm_env, error=error) + force_env%subsys => sub_force_env(primary_subsys)%force_env%subsys + CALL cp_subsys_retain(force_env%subsys, error=error) + + ENDIF + + END SUBROUTINE qmmm_update_force_mixing_env + ! ***************************************************************************** !> \param error variable to control error logging, stopping,... !> see module cp_error_handling @@ -83,18 +266,18 @@ CONTAINS !> 05.2004 created [fawzi] !> \author Fawzi Mohamed ! ***************************************************************************** - SUBROUTINE qmmm_create_force_env(force_env, root_section, para_env, globenv,& - force_env_section, subsys_section, use_motion_section, error) + SUBROUTINE qmmm_create_force_env_low(force_env, root_section, para_env, globenv,& + force_env_section, qmmm_section, subsys_section, use_motion_section, error) TYPE(force_env_type), POINTER :: force_env TYPE(section_vals_type), POINTER :: root_section TYPE(cp_para_env_type), POINTER :: para_env TYPE(global_environment_type), POINTER :: globenv TYPE(section_vals_type), POINTER :: force_env_section, & - subsys_section + qmmm_section, subsys_section LOGICAL, INTENT(IN) :: use_motion_section TYPE(cp_error_type), INTENT(inout) :: error - CHARACTER(len=*), PARAMETER :: routineN = 'qmmm_create_force_env', & + CHARACTER(len=*), PARAMETER :: routineN = 'qmmm_create_force_env_low', & routineP = moduleN//':'//routineN CHARACTER(len=default_string_length), & @@ -126,7 +309,7 @@ CONTAINS TYPE(qmmm_links_type), POINTER :: qmmm_links TYPE(qs_environment_type), POINTER :: qs_env TYPE(section_vals_type), POINTER :: print_gen, print_section, & - qmmm_periodic, qmmm_section + qmmm_periodic CALL timeset(routineN,handle) failure=.FALSE. @@ -146,7 +329,6 @@ CONTAINS subsys_section => section_vals_get_subs_vals(force_env_section,"SUBSYS",& error=error) END IF - qmmm_section => section_vals_get_subs_vals(force_env_section,"QMMM",error=error) qmmm_periodic => section_vals_get_subs_vals(qmmm_section,"PERIODIC",error=error) print_section => section_vals_get_subs_vals(qmmm_section,"PRINT",error=error) print_gen => section_vals_get_subs_vals(print_section,"PROGRAM_RUN_INFO",error=error) @@ -209,12 +391,12 @@ CPSourceFileRef,& ! First Initialize Fist... CALL section_vals_val_set(force_env_section,"METHOD",i_val=do_fist,error=error) - CALL fist_create_force_env(sub_force_env(1)%force_env, & + CALL fist_create_force_env(sub_force_env(fist_subsys)%force_env, & root_section, para_env, globenv, qmmm=.TRUE., qmmm_env=qmmm_env_mm,& force_env_section=force_env_section, subsys_section=subsys_section,& use_motion_section=use_motion_section,error=error) - CALL force_env_get(sub_force_env(1)%force_env,subsys=subsys_mm,& + CALL force_env_get(sub_force_env(fist_subsys)%force_env,subsys=subsys_mm,& cell=mm_cell,fist_env=fist_env,error=error) ! Set up QM/MM Options @@ -328,13 +510,13 @@ CPSourceFileRef,& force_env_section=force_env_section, subsys_section=subsys_section,error=error) IF (qmmm_link_imomm) CALL qmmm_link_Imomm_coord(qmmm_links, subsys_qm%particles%els,& qm_atom_index, error) - CALL quickstep_create_force_env(sub_force_env(2)%force_env,root_section,para_env,& + CALL quickstep_create_force_env(sub_force_env(qs_subsys)%force_env,root_section,para_env,& globenv, subsys=subsys_qm, cell=qm_cell_small, qmmm=.TRUE.,qmmm_env_qm=qmmm_env_qm,& qmmm_periodic=qmmm_env_qm%periodic, force_env_section=force_env_section,& subsys_section=subsys_section,use_motion_section=use_motion_section,& error=error) CALL cp_subsys_release(subsys_qm,error=error) - CALL force_env_get(sub_force_env(2)%force_env,qs_env=qs_env,error=error) + CALL force_env_get(sub_force_env(qs_subsys)%force_env,qs_env=qs_env,error=error) IF (qmmm_env_qm%periodic) THEN IF (.NOT.ASSOCIATED(super_cell)) THEN ALLOCATE (super_cell,STAT=stat) @@ -414,6 +596,6 @@ CPSourceFileRef,& END IF CALL timestop(handle) - END SUBROUTINE qmmm_create_force_env + END SUBROUTINE qmmm_create_force_env_low END MODULE qmmm_main diff --git a/src/qmmm_types.F b/src/qmmm_types.F index a2dead8910..1a0d7073fc 100644 --- a/src/qmmm_types.F +++ b/src/qmmm_types.F @@ -35,8 +35,17 @@ MODULE qmmm_types LOGICAL, PRIVATE, PARAMETER :: debug_this_module=.TRUE. CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qmmm_types' INTEGER, SAVE, PRIVATE :: last_qmmm_env_id_nr=0 + INTEGER, PARAMETER, PUBLIC :: primary_subsys = 1 INTEGER, PARAMETER, PUBLIC :: fist_subsys = 1,& qs_subsys = 2 + INTEGER, PARAMETER, PUBLIC :: force_mixing_extended_subsys = 1,& + force_mixing_core_subsys = 2 + INTEGER, PARAMETER, PUBLIC :: force_mixing_label_none = -1,& + force_mixing_label_QM_core_list = 10,& + force_mixing_label_QM_core = 9,& + force_mixing_label_QM_extended = 8,& + force_mixing_label_buffer = 7,& + force_mixing_label_termination = 6 PUBLIC :: qmmm_env_qm_type, gridlevel_info_type, qmmm_pot_type, qmmm_pot_p_type PUBLIC :: qmmm_env_qm_retain, qmmm_env_qm_release, qmmm_env_qm_create diff --git a/src/qmmm_util.F b/src/qmmm_util.F index 24ffa2f0e3..50a6caa505 100644 --- a/src/qmmm_util.F +++ b/src/qmmm_util.F @@ -9,13 +9,15 @@ !> \author Teodoro Laino ! ***************************************************************************** MODULE qmmm_util - USE cell_types, ONLY: cell_type + USE cell_types, ONLY: cell_copy,& + cell_type USE cp_subsys_types, ONLY: cp_subsys_type USE f77_blas USE force_env_types, ONLY: force_env_get,& force_env_type,& use_qmmm - USE input_constants, ONLY: do_qmmm_wall_quadratic,& + USE input_constants, ONLY: do_qmmm_wall_none,& + do_qmmm_wall_quadratic,& do_qmmm_wall_reflective USE input_section_types, ONLY: section_vals_get,& section_vals_get_subs_vals,& @@ -23,8 +25,12 @@ MODULE qmmm_util section_vals_val_get USE kinds, ONLY: dp USE particle_types, ONLY: particle_type,& + write_fist_particle_coordinates,& write_qs_particle_coordinates USE qmmm_types, ONLY: fist_subsys,& + force_mixing_core_subsys,& + force_mixing_extended_subsys,& + primary_subsys,& qs_subsys USE qs_energy_types, ONLY: qs_energy_type USE qs_environment_types, ONLY: get_qs_env @@ -33,11 +39,12 @@ MODULE qmmm_util IMPLICIT NONE PRIVATE - LOGICAL, PRIVATE, PARAMETER :: debug_this_module=.TRUE. + LOGICAL, PRIVATE, PARAMETER :: debug_this_module=.FALSE. CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qmmm_util' PUBLIC :: apply_qmmm_walls_reflective,& apply_qmmm_walls,& apply_qmmm_translate,& + apply_qmmm_translate_low,& spherical_cutoff_factor CONTAINS @@ -119,6 +126,8 @@ CONTAINS CALL section_vals_val_get(walls_section,"TYPE",i_val=iwall_type,error=error) skin(:) = list(:) ELSE + ![NB] + iwall_type=do_qmmm_wall_reflective skin(:) = 0.0_dp END IF CPPrecondition(ASSOCIATED(force_env),cp_failure_level,routineP,error,failure) @@ -126,9 +135,16 @@ CONTAINS CPPrecondition(ASSOCIATED(force_env%qmmm_env),cp_failure_level,routineP,error,failure) CPPrecondition(force_env%qmmm_env%ref_count>0,cp_failure_level,routineP,error,failure) - CALL force_env_get(force_env%sub_force_env(fist_subsys)%force_env,& + CALL cp_assert(iwall_type == do_qmmm_wall_none .OR. SIZE(force_env%sub_force_env) == 1, & + cp_warning_level,cp_assertion_failed,routineP,& + "Reflective walls for QM/MM are not implemented (or useful), when "//& + "top level force_env has two or more sub_force_env allocated, which only happens when "//& + "force mixing is active. Skipping!"//CPSourceFileRef) + IF (SIZE(force_env%sub_force_env) /= 1) RETURN + + CALL force_env_get(force_env%sub_force_env(primary_subsys)%force_env%sub_force_env(fist_subsys)%force_env,& cell=mm_cell,subsys=subsys_mm,error=error) - CALL force_env_get(force_env%sub_force_env(qs_subsys)%force_env,& + CALL force_env_get(force_env%sub_force_env(primary_subsys)%force_env%sub_force_env(qs_subsys)%force_env,& cell=qm_cell,subsys=subsys_qm,error=error) qm_atom_index => force_env%qmmm_env%qm_atom_index CPPrecondition(ASSOCIATED(qm_atom_index),cp_failure_level,routineP,error,failure) @@ -228,10 +244,19 @@ CONTAINS CPPrecondition(force_env%ref_count>0,cp_failure_level,routineP,error,failure) CPPrecondition(ASSOCIATED(force_env%qmmm_env),cp_failure_level,routineP,error,failure) CPPrecondition(force_env%qmmm_env%ref_count>0,cp_failure_level,routineP,error,failure) - CALL force_env_get(force_env%sub_force_env(fist_subsys)%force_env,& + + CALL cp_assert(SIZE(force_env%sub_force_env) == 1, & + cp_warning_level,cp_assertion_failed,routineP,& + "Quadratic walls for QM/MM are not implemented (or useful), when "//& + "top level force_env has two or more sub_force_env allocated, which only happens when "//& + "force mixing is active. Skipping!"//CPSourceFileRef) + IF (SIZE(force_env%sub_force_env) /= 1) RETURN + + CALL force_env_get(force_env%sub_force_env(primary_subsys)%force_env%sub_force_env(fist_subsys)%force_env,& cell=mm_cell,subsys=subsys_mm,error=error) - CALL force_env_get(force_env%sub_force_env(qs_subsys)%force_env,& + CALL force_env_get(force_env%sub_force_env(primary_subsys)%force_env%sub_force_env(qs_subsys)%force_env,& cell=qm_cell,subsys=subsys_qm,error=error) + qm_atom_index => force_env%qmmm_env%qm_atom_index CPPrecondition(ASSOCIATED(qm_atom_index),cp_failure_level,routineP,error,failure) @@ -294,11 +319,11 @@ CONTAINS ENDIF ENDDO - CALL get_qs_env(qs_env=force_env%sub_force_env(qs_subsys)%& + CALL get_qs_env(qs_env=force_env%sub_force_env(primary_subsys)%force_env%sub_force_env(qs_subsys)%& force_env%qs_env, energy=energy,error=error) energy%total = energy%total + wallenergy - END SUBROUTINE + END SUBROUTINE apply_qmmm_walls_quadratic ! ***************************************************************************** !> \brief Apply translation to the full system in order to center the QM @@ -316,9 +341,95 @@ CONTAINS CHARACTER(len=*), PARAMETER :: routineN = 'apply_qmmm_translate', & routineP = moduleN//':'//routineN + INTEGER :: ip + INTEGER, DIMENSION(:), POINTER :: qm_atom_index + LOGICAL :: failure + REAL(dp), POINTER :: charges(:) + TYPE(cell_type), POINTER :: cell_core, cell_extended + TYPE(cp_subsys_type), POINTER :: subsys_core, subsys_extended, & + subsys_fist, subsys_qs + TYPE(particle_type), DIMENSION(:), & + POINTER :: particles_core, & + particles_extended, & + particles_mm, particles_qm + TYPE(section_vals_type), POINTER :: subsys_section + + SELECT CASE(force_env%in_use) + + CASE(use_qmmm) + IF (SIZE(force_env%sub_force_env) == 1) THEN + CALL apply_qmmm_translate_low(force_env%sub_force_env(1)%force_env,error) + ELSE IF (SIZE(force_env%sub_force_env) == 2) THEN + ! want to center extended, and make core consistent with that + CALL apply_qmmm_translate_low(force_env%sub_force_env(force_mixing_extended_subsys)%force_env,error) + + ! translate core fist particles + CALL force_env_get(force_env%sub_force_env(force_mixing_extended_subsys)%force_env,& + subsys=subsys_extended,cell=cell_extended,error=error) + CALL force_env_get(force_env%sub_force_env(force_mixing_core_subsys)%force_env,& + subsys=subsys_core,cell=cell_core,error=error) + particles_extended => subsys_extended%particles%els + particles_core => subsys_core%particles%els + DO ip=1,SIZE(particles_extended) + particles_core(ip)%r = particles_extended(ip)%r + END DO + CALL cell_copy(cell_extended, cell_core, error) + + ! also update core QM particles (like in apply_qmmm_translate_low) + CALL force_env_get(force_env%sub_force_env(force_mixing_core_subsys)%force_env%sub_force_env(fist_subsys)%force_env,& + subsys=subsys_fist,error=error) + CALL force_env_get(force_env%sub_force_env(force_mixing_core_subsys)%force_env%sub_force_env(qs_subsys)%force_env,& + subsys=subsys_qs,error=error) + qm_atom_index => force_env%sub_force_env(force_mixing_core_subsys)%force_env%qmmm_env%qm_atom_index + particles_qm => subsys_qs%particles%els + particles_mm => subsys_fist%particles%els + CPPrecondition(ASSOCIATED(qm_atom_index),cp_failure_level,routineP,error,failure) + DO ip=1,SIZE(qm_atom_index) + particles_qm(ip)%r=particles_mm(qm_atom_index(ip))%r + END DO + +IF (debug_this_module) THEN +CALL force_env_get(force_env%sub_force_env(force_mixing_core_subsys)%force_env%sub_force_env(qs_subsys)%force_env,& + subsys=subsys_core,cell=cell_core,error=error) +particles_core => subsys_core%particles%els +subsys_section => section_vals_get_subs_vals(force_env%force_env_section, & + "SUBSYS",error=error) +CALL write_qs_particle_coordinates(particles_core,subsys_section,"QM/MM core calc first QM, then MM (0 charges)",error) +CALL force_env_get(force_env%sub_force_env(force_mixing_core_subsys)%force_env%sub_force_env(fist_subsys)%force_env,& + subsys=subsys_core,cell=cell_core,error=error) +particles_core => subsys_core%particles%els +ALLOCATE(charges(SIZE(particles_core))) +charges = 0.0_dp +CALL write_fist_particle_coordinates(particles_core,subsys_section,charges,error) +DEALLOCATE(charges) +ENDIF + + ELSE + CPPrecondition(.FALSE.,cp_failure_level,routineP,error,failure) + ENDIF + END SELECT + END SUBROUTINE apply_qmmm_translate + +! ***************************************************************************** +!> \brief Apply translation to the full system in order to center the QM +!> system into the QM box +!> \param error variable to control error logging, stopping,... +!> see module cp_error_handling +!> \par History +!> 08.2007 created [tlaino] - Zurich University +!> \author Teodoro Laino +! ***************************************************************************** + SUBROUTINE apply_qmmm_translate_low(force_env,error) + TYPE(force_env_type), POINTER :: force_env + TYPE(cp_error_type), INTENT(inout) :: error + + CHARACTER(len=*), PARAMETER :: routineN = 'apply_qmmm_translate_low', & + routineP = moduleN//':'//routineN + INTEGER :: ip, unit_nr INTEGER, DIMENSION(:), POINTER :: qm_atom_index LOGICAL :: failure + REAL(DP), POINTER :: charges(:) REAL(KIND=dp), DIMENSION(3) :: max_coord, min_coord, transl_v TYPE(cell_type), POINTER :: mm_cell, qm_cell TYPE(cp_logger_type), POINTER :: logger @@ -383,12 +494,18 @@ CONTAINS DO ip=1,SIZE(qm_atom_index) particles_qm(ip)%r=particles_mm(qm_atom_index(ip))%r END DO + subsys_section => section_vals_get_subs_vals(force_env%force_env_section, & "SUBSYS",error=error) - CALL write_qs_particle_coordinates(particles_qm,subsys_section,"QM/MM",error) + CALL write_qs_particle_coordinates(particles_qm,subsys_section,"QM/MM first QM, then MM (0 charges)",error) + ALLOCATE(charges(SIZE(particles_mm))) + charges = 0.0_dp + CALL write_fist_particle_coordinates(particles_mm,subsys_section,charges,error) + DEALLOCATE(charges) + END SELECT - END SUBROUTINE apply_qmmm_translate + END SUBROUTINE apply_qmmm_translate_low ! ***************************************************************************** !> \brief Computes a spherical cutoff factor for the QMMM interactions diff --git a/tests/QMMM/QS/regtest-1/H2O-qmmm-none-force-mixing-1.inp b/tests/QMMM/QS/regtest-1/H2O-qmmm-none-force-mixing-1.inp new file mode 100644 index 0000000000..f1153e68fd --- /dev/null +++ b/tests/QMMM/QS/regtest-1/H2O-qmmm-none-force-mixing-1.inp @@ -0,0 +1,122 @@ +&FORCE_EVAL + METHOD QMMM + &DFT + BASIS_SET_FILE_NAME ../BASIS_SET + POTENTIAL_FILE_NAME ../POTENTIAL + &MGRID + COMMENSURATE + CUTOFF 50 + &END MGRID + &QS + &END QS + &SCF + SCF_GUESS atomic + &END SCF + &XC + &XC_FUNCTIONAL pade + &END XC_FUNCTIONAL + &END XC + &END DFT + &MM + &FORCEFIELD + parm_file_name ../sample_pot/5H2O.pot + parmtype CHM + &spline + rcut_nb 5.2917720830 + &end + &END FORCEFIELD + &POISSON + &EWALD + EWALD_TYPE spme + ALPHA .44 + GMAX 64 + O_SPLINE 6 + &END EWALD + &END POISSON + &END MM + &QMMM + # NOCENTER0 + NOCENTER + MM_POTENTIAL_FILE_NAME ../MM_POTENTIAL + &CELL + ABC 17.0 7.0 7.0 + &END CELL + ECOUPL GAUSS + USE_GEEP_LIB 6 + NOCOMPATIBILITY + + # core list + &QM_KIND H + MM_INDEX 2 3 + &END QM_KIND + &QM_KIND O + MM_INDEX 1 + &END QM_KIND + + # Noam's new section + &FORCE_MIXING + R_CORE 2.5 3.5 + R_QM 2.5 3.5 + R_BUF 2.5 3.5 + QM_KIND_ELEMENT_MAPPING H H + QM_KIND_ELEMENT_MAPPING O O +# &RESTART_INFO +# INDICES 1 2 3 4 5 6 +# LABELS 10 10 10 8 8 8 +# &END RESTART_INFO + &END FORCE_MIXING + # + # QM_KINDS + # + &END QMMM + &SUBSYS + &CELL + ABC 20.00 20.00 20.00 + &END CELL + &KIND H + BASIS_SET DZVP-GTH-BLYP + POTENTIAL GTH-BLYP-q1 + &END KIND + &KIND O + BASIS_SET DZVP-GTH-BLYP + POTENTIAL GTH-BLYP-q6 + &END KIND + &TOPOLOGY + CONN_FILE_NAME ../sample_psf/5H2O.psf + CONNECTIVITY UPSF + COORD_FILE_NAME ../sample_pdb/5H2O.pdb + COORDINATE pdb + &END TOPOLOGY + &PRINT +# &ATOMIC_COORDINATES +# &END ATOMIC_COORDINATES +# &CELL +# &END CELL + &END PRINT + &END SUBSYS + &PRINT + &FORCES + &END FORCES + &END PRINT +&END FORCE_EVAL +&GLOBAL + PROJECT H2O-qmmm-none-force_mixing-1 + PRINT_LEVEL MEDIUM + RUN_TYPE MD +&END GLOBAL +&MOTION + &MD + STEPS 8 + ENSEMBLE REFTRAJ + &REFTRAJ + TRAJ_FILE_NAME ../sample_xyz/5H2O_adaptive.xyz + EVAL_ENERGY_FORCES + &END REFTRAJ + &END MD + + &PRINT + &RESTART + &END RESTART + &END PRINT + +&END MOTION diff --git a/tests/QMMM/QS/regtest-1/TEST_FILES b/tests/QMMM/QS/regtest-1/TEST_FILES index 29c47d1221..b82ad43cce 100644 --- a/tests/QMMM/QS/regtest-1/TEST_FILES +++ b/tests/QMMM/QS/regtest-1/TEST_FILES @@ -25,3 +25,4 @@ H2O-qmmm-gauss-16.inp 7 H2O-qmmm-gauss-17.inp 7 H2O-qmmm-gauss-18.inp 7 H2O-qmmm-gauss-19.inp 7 +H2O-qmmm-none-force-mixing-1.inp 7 diff --git a/tests/QMMM/QS/sample_pdb/5H2O.pdb b/tests/QMMM/QS/sample_pdb/5H2O.pdb new file mode 100644 index 0000000000..04d5bd7000 --- /dev/null +++ b/tests/QMMM/QS/sample_pdb/5H2O.pdb @@ -0,0 +1,20 @@ +ATOM 1 O WAT 1 4.000 0.000 0.000 1.00 0.00 O +ATOM 2 H1 WAT 1 4.000 0.000 0.947 1.00 0.00 H +ATOM 3 H2 WAT 1 4.000 0.896 -0.317 1.00 0.00 H +TER +ATOM 4 O WAT 2 6.000 0.000 0.000 1.00 0.00 O +ATOM 5 H1 WAT 2 6.000 0.000 0.947 1.00 0.00 H +ATOM 6 H2 WAT 2 6.000 0.896 -0.317 1.00 0.00 H +TER +ATOM 7 O WAT 3 8.000 0.000 0.000 1.00 0.00 O +ATOM 8 H1 WAT 3 8.000 0.000 0.947 1.00 0.00 H +ATOM 9 H2 WAT 3 8.000 0.896 -0.317 1.00 0.00 H +TER +ATOM 10 O WAT 4 10.000 0.000 0.000 1.00 0.00 O +ATOM 11 H1 WAT 4 10.000 0.000 0.947 1.00 0.00 H +ATOM 12 H2 WAT 4 10.000 0.896 -0.317 1.00 0.00 H +TER +ATOM 13 O WAT 5 12.000 0.000 0.000 1.00 0.00 O +ATOM 14 H1 WAT 5 12.000 0.000 0.947 1.00 0.00 H +ATOM 15 H2 WAT 5 12.000 0.896 -0.317 1.00 0.00 H +END diff --git a/tests/QMMM/QS/sample_pot/5H2O.pot b/tests/QMMM/QS/sample_pot/5H2O.pot new file mode 100644 index 0000000000..39fee2cdb9 --- /dev/null +++ b/tests/QMMM/QS/sample_pot/5H2O.pot @@ -0,0 +1,76 @@ +*>>>>>>> AMBER FF Converted into CHARMM FF style <<<<<<< +*>>>>>>> Generated on :: 20041222 165204.289 +0100 by :: teo <<<<<<< +*>>>>>>> Teddy <<<<<<< +*>>>>>>> Leap Title :: <<<<<<< +*>>>>>>> Send all comments related to the FFs conversion to <<<<<<< +*>>>>>>> teodoro.laino@gmail.com <<<<<<< + +BONDS +! +!V(bond) = Kb(b - b0)**2 +! +!Kb: kcal/mole/A**2 +!b0: A +! +!atom type Kb b0 +! +HW OW 553.000000000 0.957200000 +HW HW 553.000000000 1.513600000 + +ANGLES +! +!V(angle) = Ktheta(Theta - Theta0)**2 +! +!V(Urey-Bradley) = Kub(S - S0)**2 +! +!Ktheta: kcal/mole/rad**2 +!Theta0: degrees +!Kub: kcal/mole/A**2 (Urey-Bradley) +!S0: A +! +!atom types Ktheta Theta0 Kub S0 +! + +DIHEDRALS +! +!V(dihedral) = Kchi(1 + cos(n(chi) - delta)) +! +!Kchi: kcal/mole +!n: multiplicity +!delta: degrees +! +!atom types Kchi n delta +! + +IMPROPER +! +!V(improper) = Kpsi(psi - psi0)**2 +! +!Kpsi: kcal/mole/rad**2 +!psi0: degrees +!note that the second column of numbers (0) is ignored +! +!atom types Kpsi psi0 +! + +NONBONDED +! +!V(Lennard-Jones) = Eps,i,j[(Rmin,i,j/ri,j)**12 - 2(Rmin,i,j/ri,j)**6] +! +!epsilon: kcal/mole, Eps,i,j = sqrt(eps,i * eps,j) +!Rmin/2: A, Rmin,i,j = Rmin/2,i + Rmin/2,j +! +!atom ignored epsilon Rmin/2 ignored eps,1-4 Rmin/2,1-4 +! +OW 0.000000000 0.152000000 1.768300000 +HW 0.000000000 0.000000000 0.000000000 + +END + +! +! This Section can be cutted & pasted into the Fist input file.. +! +CHARGES +OW -0.834000000 +HW 0.417000000 +END CHARGES diff --git a/tests/QMMM/QS/sample_psf/5H2O.psf b/tests/QMMM/QS/sample_psf/5H2O.psf new file mode 100644 index 0000000000..0583e8bc0e --- /dev/null +++ b/tests/QMMM/QS/sample_psf/5H2O.psf @@ -0,0 +1,45 @@ +PSF + + 1 !NTITLE + Conversion from AMBER PARMTOP :: + + 15 !NATOM + 1 MOL1 1 WAT O OW -0.83400 16.00000 0 + 2 MOL1 1 WAT H1 HW 0.41700 1.00800 0 + 3 MOL1 1 WAT H2 HW 0.41700 1.00800 0 + 4 MOL2 2 WAT O OW -0.83400 16.00000 0 + 5 MOL2 2 WAT H1 HW 0.41700 1.00800 0 + 6 MOL2 2 WAT H2 HW 0.41700 1.00800 0 + 7 MOL3 3 WAT O OW -0.83400 16.00000 0 + 8 MOL3 3 WAT H1 HW 0.41700 1.00800 0 + 9 MOL3 3 WAT H2 HW 0.41700 1.00800 0 + 10 MOL4 4 WAT O OW -0.83400 16.00000 0 + 11 MOL4 4 WAT H1 HW 0.41700 1.00800 0 + 12 MOL4 4 WAT H2 HW 0.41700 1.00800 0 + 13 MOL5 5 WAT O OW -0.83400 16.00000 0 + 14 MOL5 5 WAT H1 HW 0.41700 1.00800 0 + 15 MOL5 5 WAT H2 HW 0.41700 1.00800 0 + + 15 !NBOND + 3 1 2 1 2 3 6 4 + 5 4 5 6 9 7 8 7 + 8 9 10 11 10 12 11 12 + 13 14 13 15 14 15 + + 0 !NTHETA + + + 0 !NPHI + + + 0 !NIMPHI + + + 0 !NDON + + 0 !NACC + + 0 !NNB + + 0 !NGRP + diff --git a/tests/QMMM/QS/sample_xyz/5H2O_adaptive.xyz b/tests/QMMM/QS/sample_xyz/5H2O_adaptive.xyz new file mode 100644 index 0000000000..241f7de546 --- /dev/null +++ b/tests/QMMM/QS/sample_xyz/5H2O_adaptive.xyz @@ -0,0 +1,136 @@ +15 +all intermolecular distances r = 2.0 <= r_inner = 2.5, should be 6/12 +O 5.000 0.000 0.000 +H 5.000 0.000 0.947 +H 5.000 0.896 -0.317 +O 7.000 0.000 0.000 +H 7.000 0.000 0.947 +H 7.000 0.896 -0.317 +O 9.000 0.000 0.000 +H 9.000 0.000 0.947 +H 9.000 0.896 -0.317 +O 11.000 0.000 0.000 +H 11.000 0.000 0.947 +H 11.000 0.896 -0.317 +O 13.000 0.000 0.000 +H 13.000 0.000 0.947 +H 13.000 0.896 -0.317 +15 +move molecule #4 to r_inner = 2.5 < r = 3.0 <= r_outer = 3.5, still 6/12 +O 5.100 0.000 0.000 +H 5.100 0.000 0.947 +H 5.100 0.896 -0.317 +O 7.100 0.000 0.000 +H 7.100 0.000 0.947 +H 7.100 0.896 -0.317 +O 9.100 0.000 0.000 +H 9.100 0.000 0.947 +H 9.100 0.896 -0.317 +O 12.100 0.000 0.000 +H 12.100 0.000 0.947 +H 12.100 0.896 -0.317 +O 14.100 0.000 0.000 +H 14.100 0.000 0.947 +H 14.100 0.896 -0.317 +15 +move 1 atom in molecule #4 to r = 3.55 > r_outer = 3.5, rest are at r = 3.45 <= r_outer, still 6/12 +O 5.100 0.000 0.000 +H 5.100 0.000 0.947 +H 5.100 0.896 -0.317 +O 7.100 0.000 0.000 +H 7.100 0.000 0.947 +H 7.100 0.896 -0.317 +O 9.100 0.000 0.000 +H 9.100 0.000 0.947 +H 9.100 0.896 -0.317 +O 12.650 0.000 0.000 +H 12.550 0.000 0.947 +H 12.550 0.896 -0.317 +O 14.650 0.000 0.000 +H 14.650 0.000 0.947 +H 14.650 0.896 -0.317 +15 +move rest of molecule #4 to r = 4.1 > r_outer = 3.5, now 6/9 +O 5.100 0.000 0.000 +H 5.100 0.000 0.947 +H 5.100 0.896 -0.317 +O 7.100 0.000 0.000 +H 7.100 0.000 0.947 +H 7.100 0.896 -0.317 +O 9.100 0.000 0.000 +H 9.100 0.000 0.947 +H 9.100 0.896 -0.317 +O 13.200 0.000 0.000 +H 13.200 0.000 0.947 +H 13.200 0.896 -0.317 +O 15.200 0.000 0.000 +H 15.200 0.000 0.947 +H 15.200 0.896 -0.317 +15 +move one atom of molecule #4 to r_inner = 2.5 < r = 3.45 < r_outer = 3.5, rest at r = 3.55 <= r_outer = 3.5 still 6/9 +O 5.100 0.000 0.000 +H 5.100 0.000 0.947 +H 5.100 0.896 -0.317 +O 7.100 0.000 0.000 +H 7.100 0.000 0.947 +H 7.100 0.896 -0.317 +O 9.100 0.000 0.000 +H 9.100 0.000 0.947 +H 9.100 0.896 -0.317 +O 12.550 0.000 0.000 +H 12.650 0.000 0.947 +H 12.650 0.896 -0.317 +O 14.650 0.000 0.000 +H 14.650 0.000 0.947 +H 14.650 0.896 -0.317 +15 +move rest of molecule #4 to r_inner = 2.5 < r = 3.1 < r_outer = 3.5, still 6/9 +O 5.100 0.000 0.000 +H 5.100 0.000 0.947 +H 5.100 0.896 -0.317 +O 7.100 0.000 0.000 +H 7.100 0.000 0.947 +H 7.100 0.896 -0.317 +O 9.100 0.000 0.000 +H 9.100 0.000 0.947 +H 9.100 0.896 -0.317 +O 12.200 0.000 0.000 +H 12.200 0.000 0.947 +H 12.200 0.896 -0.317 +O 14.200 0.000 0.000 +H 14.200 0.000 0.947 +H 14.200 0.896 -0.317 +15 +move one atom of molecule #4 to r = 2.45 < r_inner = 2.5, rest at r_inner = 2.5 < r = 2.55 <= r_outer = 3.5, back to 6/12 +O 5.100 0.000 0.000 +H 5.100 0.000 0.947 +H 5.100 0.896 -0.317 +O 7.100 0.000 0.000 +H 7.100 0.000 0.947 +H 7.100 0.896 -0.317 +O 9.100 0.000 0.000 +H 9.100 0.000 0.947 +H 9.100 0.896 -0.317 +O 11.550 0.000 0.000 +H 11.650 0.000 0.947 +H 11.650 0.896 -0.317 +O 13.650 0.000 0.000 +H 13.650 0.000 0.947 +H 13.650 0.896 -0.317 +15 +move rest of molecule #4 to r = 2.1 < r_inner = 2.5, still 6/12 +O 5.100 0.000 0.000 +H 5.100 0.000 0.947 +H 5.100 0.896 -0.317 +O 7.100 0.000 0.000 +H 7.100 0.000 0.947 +H 7.100 0.896 -0.317 +O 9.100 0.000 0.000 +H 9.100 0.000 0.947 +H 9.100 0.896 -0.317 +O 11.200 0.000 0.000 +H 11.200 0.000 0.947 +H 11.200 0.896 -0.317 +O 13.200 0.000 0.000 +H 13.200 0.000 0.947 +H 13.200 0.896 -0.317